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ABSTRACT 


The  subset  selection  algorithm  is  extended  to 
search  for  a  best  subset  from  a  large  set  of  complex- valued 
basis  functions.  This  algorithm  is  used  to  design  digital 
finite-duration  impulse  response  (FIR)  filters  having  fewer 
coefficients  than  conventional  FIR  filters.  An  optimum 
conventional  FIR  filter  is  derived  which  has  best  uniform 
spacing  of  the  fixed  number  of  samples  which  are  to  be 
used,  and  examples  are  presented  which  show  that,  for  the 
same  number  of  coefficients,  the  complex-subsc-t-se  lection 
filter  can  give  better  results  than  the  optimum 
conventional  filter. 

The  complex  subset  selection  method  is  also  applied 
to  estimation  of  the  frequencies  of  sinusoids  in  the 
presence  of  noise.  A  windowing  technique  is  introduced  to 
increase  the  efficiency  and  accuracy  of  the  algorithm  for 
frequency  estimates.  The  results  are  compared  with 
Cramer-Rao  bounds. 
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CHAPTER  1 


INTRODUCTION 

This  is  a  report  of  work  done  on  digital 

filter  design  and  sinusoidal-tone  estimation  by  the  complex 

subset  selection  method.  In  the  complex  subset  selection 

method,  a  set  of  complex  basis  functions  are  chosen  and  a 

best  subset  is  selected  to  represent  a  desired  function. 

The  desired  functions  investigated  are  an  ideal  transfer 

function  of  a  low-pass  digital  filter,  and  the  sum  of  two 

sinusoidal  tones.  The  summed  squared-magnitude  error  is 

used  as  the  criterion  for  subset  selection. 

One  can  design  a  digital  lowpass  FIR  filter  by  use 

of  the  window  method  (1),  frequency  sampling  (2), and  other 

techniques  (3),  (4).  In  all  of  these  methods,  the  input  data 

elements  are  uniformly  spaced  in  time  with  respect  to  one 

\ 

another.  Instead,  we  search  for  the  best  H  data  elements, 
which  may  or  may  not  be  uniformly  located,  to  minimize  the 
least  squares  error  in  approximating  the  desired  transfer 
function.  Nt  compare  the  result  with  a  filter  having  the 
same  number  of  coefficients,  but  uniformly  spaced 
locations. 

He  will  also  investigate  sinusoidal-tone  parameter 
estimation,  which  is  a  very  important  subject  in 
communication  systems.  Linear  prediction  (r>)  and  spectral 
estimation  (6)  are  two  well  known  methods  for  estimating 


parameters  of  tones#  such  as  frequencies.  Both  cf  these 
methods  can  produce  biased  estimates  of  the  the  signal 
parameters  when  more  than  one  signal  is  present  (7).  V’e  use 
complex  subset  selection  to  estimate  frequencies  of  tones 
without  such  bias.  The  effect  of  additive  noise  is 
discussed#  and  the  results  are  compared  with  Cramer-Eao 
bounds  computed  by  Rife  (8).  A  windowing  technique  is 
applied  to  increase  the  accuracy  and  the  efficiency  of  the 
algorithm.  By  "windowing  technique",  we  mean  that  a  small 
set  of  basis  fuctions  is  always  chosen  over  potentially 
productive  intervals,  which  will  be  searched  for  a  best 
subset,  and  other  intervals  are  ignored.  That  is, 
potentially  productive  intervals  are  first  located,  and 
only  these  intervals  are  searched  for  the  desired 
frequencies. 

The  basis  for  complex  subset  selection  is  reviewed 
in  the  following  paragraphs: 

Regression  analysis  may  be  broadly  defined  as  the 
analysis  of  relationships  among  variables.  It  is  one  of  the 
most  widely  used  statistical  tools  because  it  provides  a 
simple  method  for  establishing  a  linear  functional 
relationship  among  variables.  The  relationship  is  expressed 
in  the  form  of  an  eguation  connecting  the  response  or 
dependent  variable  and  one  or  more  independent 

variables,  x  ,x  . . .  .  The  regression  equation  takes 

I  Z  N 


+  b  x 
N  N 


(1.1) 
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Y=  b  x  +b  x  * 
11  2  2 


where  b.b  ...... b  are  called  the  regression  coefficients 

r  8  N 

and  are  deterained  from  the  data.  However,  in  some 
applications,  one  needs  to  select  N  out  of  the  N  (H  <  N) 
(independent  variables  to  best  approximate  y  by 


y=  a  u  +a  u  +.....+a  u  (1.2) 

112  2  MM 


where  0  =  \  u  ,.....,u  \  is  a  subset  selected  from  X  = 

1  l  Mi 

•[x(  , . ,x  J.  Examining  all  (  )  =  N!/(N-M)!M!  different 

subsets  to  select  the  best  subset  involves  a  lot  of 
computation,  especially  when  N  is  large. 

Hockinq  and  Leslie  (H-L)  (10)  have  developed  an 

algorithm  for  finding  the  best  subset.  In  their  algorithm, 
they  first  test  the  variables  which  are  superior.  That  is, 
the  ones  which,  when  omitted  from  the  linear  regression 
equation,  result  in  a  large  residual  error.  If  a  best 
subset  is  obtained,  they  stop.  Otherwise,  the  variables 
with  smaller  residual  error  are  tested  until  n  best  subset 
is  obtained.  They  have  developed  a  criterion  for  examining 
only  a  fraction  of  all  possible  subsets  and  stopping  after 
the  best  subset  is  selected.  The  squared  error  criterion 
was  suggested  by  Mellow  (11)  to  be  used  in  linear 
regression  models. 

The  H-L  algorithm  was  further  investigated  by 
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LeBotte  and  Hocking  (12)  for  estimation  cf  the  linear 
regression  coefficients  to  fit  a  surface  of  the  form  in 
(1.1).  They  have  also  introduced  some  other  improvements  to 
reduce  the  amount  of  computation.  Later,  Levasseur  and 
Levis  (13)  considered  antenna  design  as  a  partial  basis 
problem.  They,  too,  applied  a  generalization  of  the  H-L 
algorithm.  The  value  of  extending  the  work  of  Levasseur  and 
Lewis  to  comple x- valued  functions  was  pointed  out  by  D.  (f. 
Tufts.  He  also  suggested  applying  the  method  to  the 
synthesis  of  digital  filters  and  to  estimation  of  signal 
parameters  in  the  presence  of  noise. 

In  the  following  chapter,  we  restate  the  methods  of 
H-L  and  Levasseur  and  Lewis  for  complex  subset  selection  to 
make  it  useful  for  the  problems  considered  in  this  thesis. 
Digital  filter  synthesis  will  be  discussed  in  chapter  3. 
And  chapter  4  is  devoted  to  estimation  of  the  frequencies 
of  sinusoidal  signals  by  the  method  of  complex  subset 


selection 


CHAPTER  2 


COMPLEX  SUBSET  SELECTION 


2. 1-  Introduction 

Our  objective  is  to  extend  an  algorithm  for  partial 
basis  selection  to  complex- valued  basis  functions.  In  this 
case,  we  have  a  large  set  of  exponentials  with  imaginary 
arguments  as  basis  functions,  the  desire  is  to  approximate 
a  given  function  by  choosing  a  best  subset  of  the  basis  set 
with  minimum  approximation  error.  The  least  squares  error 
is  used  as  the  criterion  for  subset  selection.  The 
algorithm  is  desiqned  such  that  usually  only  a  fraction  of 
all  possible  subsets  is  examined  ,  and  it  comes  to  a  halt  as 
soon  as  a  best  subset  is  selected.  Hence,  considerable 
amounts  of  time  and  computation  can  usually  be  saved  by  not 
examining  all  possible  subsets. 

This  algorithm  is  suitable  for  such  problems  as 
digital  filter  design  and  sinusiodal-tone  estimation.  This 
is  due  to  the  representation  of  the  transfer  function  of  a 
nonrecursive  digital  filter  and  the  waveform  of  multiple 
sinusoidal  signals  by  linear  combinations  of  exponential 
functions  with  imaginary  arguments. 

Some  work  has  been  done  on  partial  basis  selection 
of  real-valued  functions  by  Levasseur  and  Lewis  (13).  Their 


derivation  is  restated  in  the  next  section  for  complex 
partial  basis  selection. 

2.2-  The  Algorithm 

It  is  desired  to  approximate  a  given  function  D(x) 
by  a  best  subset  of  H  elements  chosen  from  N  linearly 
independent  complex  basis  functions  relative  to  some  error 
norm.  The  least  sguares  norm  is  discussed  because  of  its 
later  use,  but  other  norms  can  be  applied,  as  well.  , 

Let  us  define  a  set  of  N  complex  basis  functions  as 
a  vector  H,  and  a  subset  with  M  (  M  <  N)  elements  selected 
from  H  as  a  vector  U. 


T 

H  =  [h  ,h  . . . h  1  (2.2.1) 

L  1  2  NJ 

T 

U  =  Fu  ,u  . . ,u  1  '  (2.2.2) 

L  1  2 


U  C=  H 

h  (t)  =  E  XP ( Jw  t)  i=1  ,2, . ,N  (2.2.3a) 

i  i 

Or 

h  (w)  =  EXP  ( Jwt  )  i  =  1 , 2, . .  N  (2.2.3b) 

i  i 


The  subset  parameter  vector  is  defined  as 


in  which  the  components  of  G  are  distinct  and  for  some 
value  of  i 


for  some  value  of  i  when  h 


9  = 

n 


is  defined  by  (2.2.3a) 


for  some  value  of  i  when  h 


is  define  by  (2.2.3b) 


(2.  2.5) 


for  n=1,2( 


i  H  and  i~~1,2,«.«  #  N 


"T"  means  transpose,  and  J=  /7. 

The  basis  functions  can  have  two  versions  in  (2.2.3)  so 
that  the  approximation  can  be  done  either  in  the  time 
domain  or  the  frequency  domain.  One  can  replace  time  with 
position,  and  have  basis  functions  related  to  location  in 
space.  The  approximating  function  A(x)  is  of  the  form 


A  (X) 


JL  T 

=  c  u  =  C  U 


(2.2.6) 


k=1  k  k 


in  which  the  coefficient  vector  C  is  defined  by 


•  [c  . 

L  1 


The  error  of  approximation  is  defined  to  be 


E(C,G)  =  D(X)  -  i(X) 


||2=j»(x|  |d 


(xj  -  k{x)  ax 


(2.2.8) 


and  R  (x)  is  a  weight  function  which  can  be  1  or  any 
suitable  function.  The  error  E(C,G)  is  to  be  minimized  by 
choice  of  the  subset  parameter  vector  G  and  the  coefficient 
vector  C.  Choosing  G  is  equivalent  to  choosing  a  subset  of 
M  basis  functions. 

It  is  impractical  to  search  over  all  (  ^j)  subsets 

M 

to  find  the  best  one,  esecially  when  N  is  large  which  is 
the  case  in  most  situations.  Hence,  the  following  method 
is  adopted  which  is  based  on  the  Hocking  and  Leslie  (10) 
algorithm,  as  extended  by  Levasseur  and  Lewis  (13)  to  a 
partial  basis  problem. 

In  the  adopted  method  one  searches  for  the 
potential  basis  functions  which  contribute  the  most  to  the 
error  function  by  their  deletion.  That  is,  a  set  of  error 
functions,  E  ^  '  s  is  obtained  where  E^  corresponds  to  the 
error  due  to  a  best  approximation  by  using  all  basis 
functions  excluding  the  ith  one. 


=  Min.  D  (x) 


•  >  c  h 
k=1  k  k 
k^i 


(2.  2.  9) 


Then,  the  basis  functions  are  re-ordered  in  a  way  that 


their  corresponding  £. errors  are  in  ascending  order. 


E  ^  E  ^ . E  (2.2.10) 

1  2  H 


The  basis  functions  coresponding  to  the  large  E^*s  are 


usually  the  potential  basis  functions  which  are  good 
candidates  to  be  kept  at  the  begining  of  the  search.  Let  r 
be  a  set  of  r=  N  -  B  basis  functions  to  be  deleted  and  let 
ES  be  the  corresponding  approximation  error.  Also,  define  p 
to  be  initially  equal  to  r.  Now,  the  search  can  be 
implemented  in  the  following  three  steps. 

STEP-1.  Delete  the  r  basis  functions  corresponding 


to  E •  '  s  with  subscripts  less  than,  or  equal  to  r  {i  r)  , 


and  compute  the  minimum  error  ER  .  If  ER  ^  E_  ,  go  to 


P+ 


step-3.  If  E*  >  E  ,  follow  step-2. 

P+l 


STEP-2.  Increment  p  by  1,  and  compute  the  error  for 


all  (  P )  combinations  of  deleting  r  out  of  p  basis 


functions  where  deleted  basis  functions  correspond  to  E.^*s 


with  i  ^  p.  Let  ER  be  the  minimum  error  obtained  from  all 


possible  (P)  combinations.  If  ER  4  E  ,  go  to  step-3; 

r  F+l 


otherwise (ER  >  E  )  re-execute  step-2. 

F+l 


STEP-3.  If  ER  ^  E 


F-r 


a  best  subset  has  been 


achieved,  and  no  further  search  is  necessary.  Hence  the 
selected  basis  functions  are  those  which  have  been  kept  for 


best 

approximation 

with  minimum 

error  either  in 

step-1  or 

step- 

2. 

NOTE: 

Step-1  is 

to  start 

the 

search  while 

step-2  is 

executed  over  and  over  until 

a  best  subset  is  obtained. 

The  search 

algorithm  was  outlined  in 

the  above 

three 

steps,  but 

it  was 

not 

explained  why 

one  stops 

whenever  ER  X  E„  . 

S  P+l 

Now,  this 

will 

be  discussed. 

Assume  n=r  +  q  with  1  ^  q  <  H  to  be  a  general  (case. 
Let  1=  |  1 ,2,. . .  p  J  be  a  set  of  p  indices  taken  from  the  set 

|l,2, . .  N  j  with  largest  element  p.  Let  be  a  set  of  r 

indices  selected  from  I  which  result  in  minimum  error  by 
their  deletion,  and  let  I* £^1,2, ••«..,  Njbe  any  set  of  r 
indices  with  an  element  i*  >  p  +  1.  The  intention  is  to 
show  that  if 


ER  =  Kin. 

lCl 


(2.  2.  1 1) 


then , 


ER  ^ 


c  h 
i-1  i  i 
i^I  • 


E 

I  • 


(2.  2.12) 


for  any  I ' . 


The  statement  (2.2.12)  can  be  proved  as  follows: 
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ER  ^  E 

p*1 


E  ^  E  BEC AUS E  i*  >  p*  1 

p+1  i' 


Ej»  is  the  Binimua  error  obtained  by  deleting  i*th  basis 
function  plus  some  other  ones,  as  well.  Hence, 

E  <  E 

i'  I' 


and , 


ER  4  E  4  E  <  E 

p+1  i*  I' 


Thus,  at  any  point,  if  the  statement  ER  ^  ^  is  satisfied, 

no  further  search  would  be  necessary.  However,  if  EE  .<  E,  , 

r, 

is  never  satisfied,  all  the  possible  subs'- ts  will  be 
examined  to  select  the  best  cne. 

To  show  the  power  of  this  algorithm,  some  numerical 
examples  will  be  presented  in  the  following  chapters.  Also, 
a  windowing  technigue  will  be  applied  to  improve  the 
efficiency  of  the  algorithm. 


CHAPTER  3 


FILTER  SYNTHESIS 


3. 1-  Introduction 

A  finite-duration  impulse  response  (FIR)  filter  can 
be  represented  with  a  set  of  coefficients  and  delays  as 
depicted  in  Figure  3.1.  The  frequency  response  of  such  a 
filter  is  of  the  form  shown  by  equation  (3.1.1)  .  The 
object  is  to  design  an  FIR  filter  to  have  fewer 
coefficients  than  equivalent  conventional  FIR  filters, 
while  meeting  the  same  performance  objectives. 


H  (w) 


(nT) 


Figure  3.  1 
EXP (JwnT) 


(3.1.1) 


Instead  of  processing  all  the  input  data  over  a  given 
record  length  or  a  subset  of  uniformly  spaced  data 
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elements,  we  choose  the  best  M  data  elements  for  achieving 
the  desired  performance.  The  fewer  the  coefficients,  the 
fewer  multiplications  will  be  performed,  since  each 
coefficient  corresponds  to  a  multiplication  directly.  The 
complex  subset  selection  method  is  used  to  achieve  the  goal 
of  selecting  the  best  subset  of  data  elements,  from  a  given 
interval  of  input  data,  to  be  processed  by  an  FIR  filter. 


3.2-  Formulation  Of  Design 

The  desired  FIR  filter  is  a  lowpass  filter (LPF) 
havinq  a  passband,  a  stopband,  and  a  'don't  care'  (3) 
region  as  depicted  in  figure  3.2.  The  desired  frequency 
transfer  function,  D  (w)  ,  is  real  valued  and  symmetric  about 
frequency  zero.  Assurainq  that  sampled  data  is  to  be 
processed  with  a  memory  nc  greater  than  2qT  seconds,  or 
(2q+  1)  samples,  we  can  approximate  D(v)  using  all  the 
available  data  by  a  FIR  filter  with  transfer  function  A  (w) 
given  in  formula  (3.2.1). 


i  i 


i 

TsiT-p  d 


t  i— 

p  si 


Figure  3.2 


SI 


oJ 
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A  00 


a(nT) 


n=-q 


EXP  (JwnT) 


(3. 2.1) 


There  is  no  loss  of  generality  by  choosing  zero  phase  for 
D (v) ,  because  one  can  easily  synthesize  a  linear  phase 
filter  by  shifting  the  impulse  response  of  this  filter 
properly.  Thus,  the  ideal  filter  is  defined  by  equation 
(3.2.2) . 


D(w) 


i  *  p 

0  .  si  $  I W  j  «  s2 

don't  care,  otherwise 


(3.2.2) 


Let  the  observation  tine  be  kept  fixed.  That  is,  q  and  T  of 
(3.2.1)  satisfy  the  formula 


2qT  =  constant.  (3.2.3) 

Using  complex  subset  selection,  D  (w)  can  be  approximated  by 
a  subset  having  fewer  coefficients  in  the  same  observation 
time  (M  <  2q+1)  shown  by  equation  (3.2.4), 

H  (w)  =  >_h(n)  EXP  («7wT  )  (3.2.4) 
n=1  n 

where  the  's  may  or  may  not  be  uniformly  spaced.  Using 
the  method  of  complex  subset  selection,  the  large  set  of 
basis  functions  of  the  form  FXP(JwnT)  are  the  potential 
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basis  functions  of  (3.2.1)  in  which 


nth  basis  function  =  EXP  (JwnT) 


(3.2.5) 


n  =  ~q,  •  •  •  # ”  1 1  0/  1  »•  •  •  #3 


N  =  2g  +  1  . 


(3.2.6) 


A  best  subset  having  M  coefficients  is  selected  from  these 
N  basis  functions  such  that  the  least  squares  error  between 
D(w)  and  H  (w)  is  minimized,  i.e.. 


(  M  2 

Min.  Error  =  K(w)  D(w)  -  Y"  h  (n)  EXP(JvT  )  dw 

J  ifr 

•|h(n),TT1j  (3.2.7) 


where 


W(v)  = 


,  w  \  <  p  and  si  <  w  ^  s2 


,  p  <  w  <  s  1 


(3.  2.  8) 


and  the  's  are  selected  out  of  the  set 

(~  ql »  »  ..  ,  •  T,  0  ,T ,  ,  •  •  ,  q  T )  . 

From  the  orthogonality  principle  (1h)  for  optimum  choice  ,f 
the  h(n)'s,  for  given  values  cf  the  T  •  s,  it  fellows  tilt 


H  (w)  (D(w)  -  2Z  h  (n)  EX  P  (J  wT  ))  .  EXP(-JwT  )  dw  =  0. 
n  =  1  n  k 


for  k  =  1,2, 


(3.2.?) 
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which  can  be  rewritten  as  a  set  of  M  linear  equations  of 
the  form 


AH  =  B 


(3.2.10) 


with 


a  a 
11  12 

A  ”  a  • . 

21 


(3.  2.  11) 


H  = 

— 


=  h(1),h(2), . *  h  (M) 


(3.2.12) 


B  =  b  ,  b  ( 
L  1  2 


(3.  2.  13) 


where 


(2/  (T  -T  )  )  (  Sin  (s2  (T  -T  ))-Sin(s1(T  -T  )) 

i  1  i  j  i  j 

♦  Sir.  (p  (T  -T  ))  )  i/j 

i  j 

(3.  2.  14) 


2  (s2-s  1  +  p) 


(2/  (T  )  )  .  Sin  (pT  )  T 


(3.  2.  15) 


T  =0 


7W- 
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H  = 


(3.  2.  16) 


For  the  detailed  derivation  of  the  equations  (3.2.9) 

through  (3.2.17),  see  appendix  A. 

Since  A  and  B  are  both  real  valued  matrices,  the 
solution  to  the  equation  (3.2.16)  results  in  real 

coefficients  for  the  filter  (i.e.  h(n)  's  are  all  real)  . 
Solving  for  H  and  substituting  in  the  error  function  of 
(3.2.7)  is  equivalent  to  simplifying  (3.2.7)  using  (3.2.9). 

The  result  is  , 

M 

Min.  Error  =  2p  -  2^h(n)  Sin(pT  )  /T  (3.2.17) 

n=1  n  n 

h  (n)  j 


in  which  ^  h  (n)  ]  satisfies  (3.2.9),  or  equivalently 
(3.2.16). 

The  results  obtained  here  will  be  compared  with  a 
filter  having  uniformly  spaced  values  of  T  with  frequency 
response  given  by  (3.3.1),  and  we  claim  that  the  complex 
subset  selection  will  always  result  in  a  better 
approximation  to  a  desired  ideal  filter,  or  at  least  an 
approximation  as  good  as  the  one  obtained  by  the 
uniformly- spaced  method.  The  uniformly-spaced  method  is  a 
search  for  only  a  constrained  set  of  possible  solutions, 
and  many  other  possibilities  are  simply  ignored.  On  the 
other  hand,  whenever  complex  subset  selecton  is  used,  all 
possible  approximations  are  considered  .  If  it  happens  that 


•  \ 


the  best  solution  is  unifoi mly  spaced,  it  will  be  included 
as  one  of  the  possible  solutions  to  be  examined. 


3.3-  Filter  With  Uniformly -spaced  Ti^:e 
Samples 

The  filter  obtained  by  the  complex  subset  selection 

will  be  compared  with  another  filter  which  uses 

uniformly-spaced  time  samples.  For  this  filter,  the 

observation  time  (2gT)  is  fixed  and  the  number  of  the 

» 

coefficients  (M)  are  the  same  as  before  .  But,  the  basis 
function  parameters  must  be  spaced  uniformly  from  an 
initial  location  *c*  as  illustrated  in  figure  3.3. 


L  A  A  & 


C 

K- — 


Figure  3.3 


The  freguency  response  of  this  filter  is  of  the  form  shown 
by  eguation  (3.3.1),  assuming  that  (1  is  an  odd  integer  of 
the  form  8=21*1 


H  (w)  =^ 


2_  h  (n)  EX?(J  (nA+c)  w) 
n=-l  u 


(3.3.1) 


•-  .uv 
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If,  for  fixed  values  of  A.  and  c,  we  minimize  the  error 
given  by 


Error  = 


r 

W(W) 

D(w)  -  Hu(w) 

dw 


(3.3.2) 


by  choice  of  the  coefficients  h(n)  ,  We  obtain,  using  the 
orthogonality  principle  (14),  the  set  of  simultaneous 
linear  equations 


A  '  H  =  B*  ( 3(.  3 . 3 ) 

in  which  A‘  and  B'  can  be  evaluated  from  (3.2.9)  by  letting 
T^=nA  +  c  .  The  coefficients  ^h^Cn)*!  are  also  real  since  A* 
and  B'  are  real-valued  matrices. 

In  the  search  for  the  best  c,  it  was  found 
numerically  that  the  best  c  is  zero  for  the  bestZf  for  the 
cases  tried  (fcur  examples  are  illustrated  in  table  3.  1 
where  the  C  was  searched  in  the  interval  0.0  to  3.0  with 
increment  .1,  and  best  A  was  searched  in  the  interval  0.  to 
2.5  with  increment  .1).  Therefore,  we  tried  to  prove  this 
analytically.  However,  after  investigating  many  different 
methods,  we  could  OEly  show  that  C=0  is  a  good  candidate 
for  best  C,  and  could  not  prove  its  uniqueness.  One  method 
is  to  differentiate  the  error  function  with  respect  to  c, 
and  set  it  to  zero.  The  error  function  is  af  rhe  form 


-I 

where  a*  (n,n  are  elements  of  A* 


Since  A*  is  a  toeplitz 


21 


matrix  (see  equations  (3.2.11)  and  (3.2.14)  ),  it  follows 
that 


a*  (n,-m)  =  a'  (-n,m) 


(3.3.8) 


and 

dE/dc  Sin  (n^lP) /nA.  (  PCos  (-mAP)  /  (- mA) 

n=1  m=-l 

-Sin  (-mAP)  /  ( (-m£>)  **2)  ).a»(n,-m) 

♦  Sin  (-n/E) /»a.  (  P.  Cos  (mAP) /"mA 
-  Sin (m^P) /  ( (mA)  **2)  =  0.  (3.3.9) 

Thus  c=0  is  a  solution  of  (3.3.5)  ,  because  all  the  terms  in 
(3.3.7)  cancel  one  another  for  c=0.  Since  we  have  not  been 
able  to  show  analytically  the  uniqueness  of  c,  we  search 
for  the  best  c  as  we  do  for  the  best  A* 

In  the  next  section,  some  examples  are  presented  in 
order  to  explain  the  differences  of  the  two  filters. 


1 


3.4-  Results  And  Examples  ; 

i 

Presented  here  are  the  results  of  the  comparison  of  j 

the  two  filters  discussed  in  sections  3.2  and  3.3.  Many  1 

1 

examples  have  been  investigated,  and  only  two  will  be  <j 

presented  here  to  show  the  differences  in  the  two  filters.  j 

In  both  examples,  21  basis  functions  are  used  with 

2 

arguments  wnT,  (n=-g,... ,0,.. .  ,q  ,  and  T=.5) ,  and  the  | 
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subset  size  is  5  (M=5) .  The  results  are  tabulated  in  tables 
3.2  and  3.3  for  which  the  following  observations  are  made. 
Of  course,  these  observations  are  drawn  from  the 
investigation  of  a  large  set  of  examples. 

a)  The  least  squares  error  is  smaller  for  the  subset 
selection  filter  than  the  uniformly  spaced  parameter 
filter. 

b)  The  magnitude  response  is  flatter  in  the  passband 
region  in  the  subset  selection  filter  case. 

c)  The  magnitude  response  rolls  off  sharply  in  the 
transition-band  for  the  subset  selection  filter. 

a)  However,  the  maximum  sidelobe  level  is  smaller  for 
the  uniform  parameter  filter. 

The  above  observations  can  be  seen  from  tables  3.2  and  3.3 
and  also  from  magnitude  response  plots  of  figure  3.h  and 
3.5.  The  frequency  response  is  symmetric  about  frequency 
1/2T,  where  T  is  the  sampling  period.  Since  the  smallest 
sampling  period  for  21  basis  functions  is  T--.5,  all 
freguency  responses  are  plotted  over  frequencies  0  through 
1  hertz  (1/2  (.5)  =  1)  . 

It  was  also  noticed  that  the  nonzero  coefficient 
locations  are  moved  toward  (or  away  from)  the  midpoint  by 


i 


increasing  (or  decreasing)  the  passband  width.  That  is,  the 
coefficient  locations  are  around  the  midpoint  for  a  wide 
passband,  and  farther  away  for  a  narrow  passband.  The  plot 


in  figure  3. 
second  tern 


6  shows  this  effect  where  the  abscissa  is  the 


in  eguation  (3 


ffect  where 
.3.2)  (  H 


2_  h  (n)  •  Sin  (nA)  p/n/)  , 

rts-P  u 


and  the  ordinate  is  A.  .  Each  curve  is  for  a  different 


passband  value  while  each  peak  represents  the  best  A  for 
that  specific  passband  value. 


TABLE  3.  la 


M=5 

P=.  15  Si  =.  2  S2=.  5 


ERROR  FOR 

ERROR  FOP 

ERROR  FOR 

c 

A=0.  5 

A  =  1. 0 

A  =  1 .5 

0.0 

. 144201 

.092526 

.  162134 

0.  1 

. 147396 

. 093267 

.  161562 

0.2 

.156949 

. 095505 

.159873 

0.3 

.172754 

.  099286 

. 157136 

0.4 

.194671 

.  104690 

. 153471 

0.5 

. 222486 

.111818 

.  149040 

0.6 

.255916 

.  120801 

.  144  04  4 

0.7 

. 294664 

.  131788 

.  138728 

0.8 

.338348 

.144936 

.  133360 

0.9 

.386559 

.160419 

. 128241 

1.0 

.438853 

.  178403 

.123694 

1.  1 

.494751 

.  199062 

. 120051 

1.2 

.553746 

. 222544 

.117660 

1.3 

.61!  326 

.  248993 

.  116866 

1.4 

.678959 

.278519 

.11801 1 

1.5 

.  7441  13 

.  311209 

. 121422 

1.6 

.810252 

.347113 

.  127416 

1.7 

. 876849 

. 386238 

.  136  27  8 

1.8 

.943399 

. 428554 

.148267 

1.9 

1. 00941 

.473985 

.  163603 

2.0 

1. 07444 

.522403 

.  182472 

Best  C 

0.  0 

0.0 

1.  3 

ERROR 

.144201 

.  092526 

.116866 
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TABLE  3.1b 


11=5 

P=.  1  S1  =  . 

2  S2=. 5 

ERBOB  POR 

ERBOB  FOR 

c 

A=0.5 

4  =  1.  0 

0.0 

.103650 

. 032061 

0.  1 

.104927 

.032794 

0.2 

.108742 

. 034986 

0.3 

.115084 

.038536 

0.4 

. 123900 

. 043731 

0.5 

.135168 

.050260 

0.6 

.148798 

. 058206 

0.7 

.164734 

.067545 

0.8 

. 162873 

. 078256 

o.y 

.203121 

.090311 

1.0 

.225369 

.  103673 

1.1 

.249488 

.118313 

1.2 

. 275352 

. 134188 

1.3 

.302826 

. 151259 

1.4 

.  331759 

.  169477 

1.5 

.361995 

. 188798 

1.6 

. 393391 

. 209166 

1.7 

.425772 

.230528 

1.8 

.458985 

.  25  28  26 

1.  9 

.492862 

.276001 

2.0 

.527239 

. 259990 

Be  st  C 

0.  0 

0.  0 

Error 

.1 03650 

.032061 

ERROR  FOR 
^  =  1.5 

.036306 
.036345 
.036466 
. 036681 
. 037008 
.037473 
. 038109 
.  038957 
.04  0061 
.  04  1470 
.043241 
.045435 
.  04811  1 
.051339 
.  055175 
.  059  697 
.064967 
.  07 1  049 
.  07801  1 
.085911 
.  09  4  80  7 

0.  0 

. 036306 
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TABLE  3.  1c 


M=5 

P=. 15  S 1= 

.225  S2  =  .  5 

ERROR  FOR 

ERROR  FOR 

ERROR  FO 

c 

=0.5 

=1.0 

^  =  1.5 

0.0 

.093154 

.068523 

.154084 

0.  1 

.096298 

. 068842 

.  153555 

0.2 

.105711 

. 069827 

.  152005 

0.  3 

.121291 

.  071544 

.  149489 

0.4 

.146899 

.074114 

.156113 

0.5 

.170340 

.077695 

.  142028 

0.  6 

.203352 

. 082489 

.137417 

0.7 

. 241646 

. 088724 

.  132498 

0.8 

.284865 

.096650 

.  127518 

0.9 

.332622 

. 106542 

.  122751 

1.0 

.384496 

.118675 

.  11  8493 

1.  1 

. 440033 

.  133334 

. 115049 

1.2 

.498750 

.  150784 

.112741 

1.3 

.560160 

. 171278 

.  111889 

1.4 

.623744 

. 195045 

.  112813 

1.  5 

.688995 

.  222276 

.115825 

1.6 

.755392 

.253120 

.129281 

1.7 

.822430 

.  287677 

.  129281 

1.8 

.689607 

.326001 

.  140256 

1.9 

.956452 

. 368084 

.  154  369 

2.0 

1. 02250 

. 41 3856 

.171806 

Best  C 

0.0 

o 

• 

o 

1.3 

Error 

.093154 

. 068523 

. 11 1B89 

27 


TABLE  3. Id 


H=5 

P=. 1  S1=. 

225  S2=.  5 

ERROR  FOR 

ERROR  FOP 

ERROR  FOR 

c 

A  =0.5 

A=1.0 

A  =  1.5 

0.0 

.062489 

.018249 

.036055 

0.  1 

.063793 

. 018802 

. 036085 

0.2 

.067678 

. 020460 

.036176 

0.3 

.074139 

. 023223 

.  036341 

0.4 

. 083124 

. 027090 

.036599 

0.5 

.094604 

. 032058 

.036975 

0.6 

. 108495 

.038125 

.037506 

0.7 

.124734 

.045289 

.038225 

0.8 

.  143224 

. 053546 

.039182 

0.9 

.263868 

.062889 

.040426 

1.0 

. 186553 

. 073311 

.042016 

1.1 

.211153 

. 004809 

.044006 

1.2 

.237535 

. 097367 

.  046464 

1.3 

.265566 

.  110980 

.049  453 

1.4 

. 395094 

.  125629 

.053040 

1.5 

.325968 

.141302 

. 057292 

1.6 

. 358029 

. 157980 

.062277 

1  7 

.391109 

. 175643 

.068061 

1.8 

. 425056 

. 194265 

.  07471 1 

1.9 

.459692 

.  213823 

.082287 

2.0 

.494855 

.  2  34  284 

.090  64  7 

Best  C 

o 

• 

o 

0.0 

0.0 

Error 

.062489 

.  018249 

.  03  6  055 

TABLE  3.2 


I 

N=  2 1  M=5 

V-15  h> 

Filter 
with  H 

coefficients 

-.2 

Uniform 

parameter 

filter 

coefficient 

nT 

locations 

n=-10 ,. . .  ,10 

0,±1.  0,12.0 

Hax.  passband 
level  (db) 

.677 

1.  3 

magnitude  (db) 
at  fp 

-2.6 

-4.  5 

magnitude  (db) 

at  H\ 

-12.9 

-10.  1 

Kax.  sidelobe 
level  (db) 

-25.  2 

-24.  3 

least  squares 
error 

.  209 8 836 E- 01 

.  9  2527E-0  1 

Tn 


Subset 
selectio  n 
filter 

Of  ±1.5, ±4. 5 
.73 

-3.5 

-17.1 

-21.0 

74  60785E-0  1 


NCY-hZ 


a 


ency  response  of 
avine  5  coeffici 


ncy  response  of  sub 
ving  5  coefficients 

f „o-.5  • 
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CHAPTEB  4 

ESTIMATION  OF  THE  FREQUENCIES 
OF  SINUSOIDS 

4.  1-  Introduction 

This  chapter  is  devoted  to  the  problem  of 
estimation  of  the  frequencies  of  sinusoidal  tones  by  the 
complex  subset  selection  method.  In  this  case,  the  function 
to  be  approximated  contains  some  tones  with  specific 
amplitudes  and  frequencies  in  the  presence  cr  absence  of 
noise.  He  consider  the  case  of  two  sinusoidal  tones,  but 
the  method  can  be  extended  to  more  tones  at  the  cost  of 
more  computation. 

The  framework  of  this  problem  will  be  developed  in 
section  4.2;  while  the  frequency  estimation  will  be 
presented  in  section  4. 3  .  The  presence  of  Gaussian  noise 
in  the  data  will  be  investigated  in  section  4.4,  and 
finally  section  4.5  will  contain  some  numerical  results. 

4.2-  Two  Tone  Approximation 

The  signal  tc  be  approximated  contains  two 
sinusiodal  tones  with  each  tone  having  amplitude  and 
frequency  (A1,F1)  and  (A2,F2)  ,  respectively.  A1  and  A2  can 


be  either  real  cr  complex 


S(t)  =  A 1  EXP  (J2j7f  t)  «•  A2  EXP  ( J2T7f  t)  (4.2.1) 

1  2 

A  data  vector  y  is  formed  from  N  samples  of  above  function 
(4.2.2)  . 


y  =  [y (0)  ,y  (1) , . y(N-i)J  (4.2.2) 

y(n)  =  S(nT) 

t 

The  complex  exponential  basis  functions  used  here  are 
h  (t)  ,h  (t)  , . . h  (t)  with  arguments  of  the  form  2flf  t  , 

I  2  L  1 

i.  e. , 

h  (t)  =  EXP(J2fff  t)  for  I- 1  ,2  ,....,  L  (4.2.3) 

I  I 

A  vector  is  formed  for  each  basis  function  with  N  samples 
taken  from  each  of  them  . 


T 

H 

I 


h  (0),h  (1),....,h  (N-1) 

II  I 


for  I=1,2,...,L 


with 


(4.2.4) 


h  (n)  =  EXP(J2Ttf  n) 

I  I 

and 

0  <  f  <  f  < .  <  f  <  1 
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Due  to  the  periodicity  of  exponentials  with  an  imaginary 
arguments,  f  should  be  less  than  one.  That  is 


EXP(J2T7fLn)  =  EXP  (J2T7(f^-1)  n) 


Having  H  ^  ‘s  for  basis  functions,  we  want  to 
approximate  the  data  vector  Y  by  the  complex  subset 
selection  method  with  minimum  least  squares  error.  In 
another  words,  we  want  to  choose  M  functions  from  the  L 
basis  functions  to  best  represent  the  data  vector,  Y  in 
order  to  minimize  the  sum  cf  the  squares  of  the  errors. 
Assume  Y  is  the  data  approximation  vector  of  the  form 


Y  =  c  u  =  UC 
k=1  k  k 


(4.2.5) 


where  U  is  a  vector  formed  from  H-column  vectors  U  's  and  C 

K 

is  a  M-dimensional  complex  coefficient  vector. 


T 

C  = 


=  ]c 

u  1 


(4.2.6) 


u  =  f o  ;u  '. . lu  1 

L  1  •  2*  *  H  J 


(4.2.7) 


T 

U  =  fu  (0)  ,u  (1)  . . U  ( N- 1 )  1 

k  L  k  k  k  J 


(4.2.3) 


u  (n)  =  EXP  (Jw  n) 
k  k 


(4.2.9) 


-  Cv 


Note  0  *s  are  chosen  from  H  *s  ,  and  v  *s  are  selected  from 

K  IK 

2  J7f^*s.  The  error  function  is 


Bin.  Error  =  E  =  j  j  Y  -  Y~ 

{c,.} 


(4. 2. 10) 


where  ;i*j  means  the  L2  norm  of  the  vector. 

To  minimize  the  error  one  has  to  find  the  optimum 
coefficient  vector  and  frequency  vector  to  be  used  in 

A 

specifying  Y. 

For  the  optimum  coefficient  vector  ,  the 


ort  oqonality  principle  indicate;  that  Y- Y  must  be 

orthogonal  to  every  Uj  with  1=1/ . ,H  .  (*  mean j  conjugate 

transpose) . 


0  (Y  -  Y)  =  0 
1 


V  * 

0  Y  =  U  Y 


U  0C  =  0  I 
1  1 


for  1=1,2. 


for  1=1,2, . .  M 


for  1=1,2, 


Let  A  =  0*0  ,  and  B  =  0* Y,THEN 


AC  =  B 


(4.  2.  11) 


(4. 2.12) 


(4.  2.  13) 


(4.2.  14) 
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where 


*  *  * 

DU  DO . D  U 

11  12  1  M 

A  =  •  «  • 

*  * 

D  0 . .  .  .  .  U  0 

Ml  MM 


with 


*  N-1 

DO  =  EXP  (J  (w  -w  )  n) 
1  k  n=0  k  1 


T  r*  *  *  * 

b  =  o  y,u  y, . ,u  y 

L  1  2  M  J 


(4.  2.  15) 


*  N- 1 

b  =  U  Y  =  y(n)  EXP  (-Jw  r.)  (4.2.16) 

1  1  n=0  1 


The  solution  of  the  linear  equation  (4.2.14)  will  give  the 
optimum  coefficient  vector  C.  Thus,  the  error  function 
becomes 


•b 

- 

^-|2 

Y II  “ 

♦ 

(Y  -  ¥>)  (Y  -  Y) 

* 

A,*  A 

I  (Y  -  Y) 

*  Y 

(I 

“  I)  - 

* 

A 

*  *  /S 

=  Y 

<Y 

'  Y)  - 

CD  (Y  -  Y) 

(4.2.17) 


The  right  side  of  the  above  equation  can  be  simplified  by 


•.*  jMtV> 


using  equation  (4.2.11) 


Y(X-/Y)  =  |y|  -  Y  'y 

2  *  2  * 

|  Y  |  -  Y  DC  =  |y  /  -  EC 


(4. 2. 18) 


2  -i2  *  N - 1  ,  . 2  M 

E  -  lY  -  B  C  =y_  (  y  (n)  .  -Y~c  y  (Q)  EXP(-Jw  n)  ) 
11  {^0  ■  '  tfn  k  k 

(«. 2.  19) 


Let  us  choose  the  frequencies  f^*s  of  the  basis 
functions  to  be  uniformly  spaced.  Hence  ,  the  functions 
will  be  of  the  form 


h  (n)  *  EXP  (J2  TTfln)  (4.2.20) 

I 

such  that  4  f = 1/L  is  the  difference  between  two  consecutive 
grid  points.  The  goodness  of  the  approximation  is  dependent 
on  L  and  N,  namely,,  the  number  of  the  basis  functions  and 
the  number  of  the  samples.  Since  Af=1/L  ,  any  increase  in  L 
will  result  in  having  denser  grid  points  which  in  turn  will 
increase  the  possibility  that  fl  and  f2  to  be  on  the  grid 
points.  ?1  and  f.2  being  on  the  grid  points  means  that  two 
°f  the  basis  functions  can  exactly  apprcxi mate  the  data 
r'  On  the  other  hand  ,  large  N  means  having  more  data 
points  to  process  which  wiU  also  resnlt  in  better 


approximation.  But  L  cannot  get  too  large,  because  the 
amount  of  the  computation  will  get  impractical  for  large  L. 
Thus  ,  the  windowing  technigue  described  below  will  be 
applied  whenever  very  dense  grid  points  are  needed. 


4.3-  Frequency  Estimation 

In  some  applications  one  has  to  estimate  the 
frequencies  of  the  tones  as  accurately  as  possible.  In 
those  cases,  we  set  the  number  of  the  basis  function^  M  in 
the  selected  subset  to  be  equal  the  number  of  the 
frequencies  contained  in  the  signal.  For  example,  when  we 
are  estimating  two  sinusioaal  tones  we  set  h  =  2.  But  a  very 
dense  set  of  grid  points  is  needed  in  order  for  the 
estimation  to  be  accurate,  and  having  a  very  dense  grid 
points  is  impractical  computationally  since  the  amount  of 
the  computation  will  increase  significantly  as  the  number 
of  the  grid  points  becomes  large.  To  qet  around  this 
deficiency,  we  introduce  the  windowing  technigue  below. 

The  windowing  technique-  is  applied  whenever  the 
grid  points  (number  of  basis  functions)  are  large,  and  the 
search  for  a  best  subset  needs  too  much  computer 
processing.  In  those  casts,  a  best  subset  is  selected  from 
a  less  dense  grid  as  a  first  step.  Ir.  second  step,  a  denser 
arid  is  chosen  around  the  selected  subset  grid  points  of 
step  one.  As  step  three,  a  new  best  subset  is  selected  from 
the  new  grid  points  of  step  two.  Then,  step  two  and  three 
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repeated  until  no  changes  are  noticed  in  the  selected 
subset  and  the  error  function.  One  should  notice  thac  in 
the  first  attempt,  Dotential  basis  functions  are  selected 
where  in  the  second  attempt  a  new  set  of  basis  functions 
are  established  which  are  densely  populated  around  those 
selected  potential  basis  functions,  and  the  best  subset  is 
selected.  Graphically,  the  windows  in  Figure  4.1  indicate 
where  the  new  densely  populated  grid  points  are  to  be  set. 


'  1 
j 

D  .1  ■ l  -5  •'i  -v  l 

'  i  •  D  •  i  1 

Figure  4. 1 


To  clarify  tie  windowing  technique,  let  us  take  a 
lock  at  the  two-tcne  siqual,  with  frequencies  f1-.2S  and 
f2=.725.  Choose  only  ten  grid  points  uniformly  spaced  iron 
zero  to  one,  for  the  first  step.  Thus,  the  basis  functions 
are 

h  (n)  =  EXP  (J2  TT f I n )  1=0,1,....  ,9 

I 

with  L=1 0  ,and  grid  points 

(0.  ,0.  1,0.  2, ......  0.9) 


of  which  2  are  selected  as  potential  functions.  Assume 
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pth=.3  and  gth=.7  are  selected.  Tn  the  next  step,  we  use 
the  following  ten  new  basis  functions  around  the  pth,  and 
qth  : 


h  (n)  =  EXP(J2TT(P*  -y-d-2)  )  n) 
I  L 


h  <n)  =  EXP  (i32'fl(q  +  4£-(l“7)  )  n) 
I  c 


1=0, 1,2,  3,  4 


1=5, 6, 7, 8,9 


with  L=1 0,  and  grid  points 


(.  2,.  25,.  3,  .3  5, .  4,  .  6,  .65,. 7  ,.75,.  8) 


Note  the  new  grid  points  are  one-half  of  the  previous  step. 
The  complex  subset  selection  algorithm  is  again  used  to 
select  two  potential  functions  from  these  new  basis 
functions.  In  this  case,  second  and  nineth  (p=.  25 , g=. 75) 
are  selected.  Again,  choose  another  ten  basip  functions 
with  grid  points  around  .25, and  .75 


h  (n)  =  EXP(J2TT(P  +  4^<1'2)  )n) 

I  * 


h  (n)  =  EXP  (J2Tr<g  +  T.-<1-7n  n) 
I  * 


with  L=10,  and  new  grid  points 


I-  0, 1,..  4 


15,6,.. 9 


(.2,.  2 25,. 25,. 275,. 3,. 67 5,. 7,. 72 5,. 75, .775)  . 


Continue  as  before  until  the  selected  frequencies  and  the 
error  do  not  change  very  much  for  a  denser  grid  points.  In 
this  case  .25  and  .725  are  the  selected  frequencies  .  The 
above  example  is  listed  in  Table  4.4. 

If  at  any  time,  two  windows  overlap,  that  is,  two 
selected  frequencies  are  close  to  one  another,  use  one  wide 
window  instead  of  two  as  depicted  below. 


— I - -  -  -  _t .. 

\  '1 

Figure  4.2 


The  program  listings  for  the  two  frequency  estimation  by 
subset  selection  and  windowing  xe-chnique  are  in  appendix  B. 


4.4-  The  Effect  Of  Additive  Noise 

The  approximation  of  a  two-tone  signal  was 
investigated  in  section  4.2.  But,  it  was  for  deterministic 
signals,  and  in  almost  all  practical  applications  the 
received  signal  is  purturbed  by  some  kind  of  noise.  In 
here,  we  investigate  the  additive  Gaussian  noise  which  is 


the  general  case  for  most  systems. 
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Recall  the  data  vector  Y  (4.2.2)  which  was  obtained 
from  a  deterministic  signal,  only.  Now  we  form  a  new  data 
vector  from  signal  vector  S  with  noise  vector  V  added  to 
it. 


Y  =  S  +  V  (4.4.1) 

T 

v  =rvC°)^vn)» . ,  V  (V- 1)  "1  (4.4.2) 


Note  that  each  vector  in  (4.4.1)  is  N- dimensiona  1.  Assume  V 

i 

consists  of  complex,  independent,  zero-mean  Gaussian  noise 
components  with  known  variance. 

By  a  close  look  at  linear  equation  (4.2.1)  ,  one 
will  realize  that  the  ncise  will  affect  the  B  vector  , 
only.  That  is,  each  rcw  of  B  is  changed  by  a  amount  due  to 
the  inner  product  of  the  noise  vector  and  the  corresponding 
U  vector  (where  subscript  corresponds  to  I the  row.) 

*  *  * 

B=UY=US*UV  (4.4.3) 

Therefore,  the  same  procedure  as  before  is  used  to  find  the 
coefficient  vector  ,  selected  basis  fuction  parameters, an d 
the  error  function  with  new  data  vector. 

The  effects  of  Gaussian  ncise  in  tone  parameter 
estimation  has  been  well  investigated  by  David  c.  Rife  (8). 
He  has  used  a  multiparameter  generalization  of  the 


Cramer-Rao  bound  (C-R  bound)  to  establish  lower  bounds  for 
variances.  He  has  also  shown  that  the  single-tone  bounds 
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are  the  minimum  values  attained  by  two-tone  bounds,  and 
they  are  attained  wher  the  two  tones  have  wide  frequency 
separation.  Thus,  one  can  use  the  single-tone  bounds 
(4.4.4)  for  two-tone  bounds  whenever  two  frequencies  are 
far  apart. 

4  2  2  2  2 

var  (w  )  >  (12  6  )/(b  T  N(N  -1))  (4.4.4) 

i  i 

where 

t 

w  frequency  estimate  in  radians 

1 

2 

6  noise  variance 

b  tone  amplitude 

i 

T  time  between  two  samples 

N  number  of  samples 

In  the  next  section,  this  C-R  bound  is  compared 
with  the  measured  accuracy  .in  frequency  estimation  which  is 
obtained  by  applyinq  the  complex  subset  selection  algorithm 
to  simulated  data. 

4.5  Results  and  Examples 

Presented  here  are  numerical  results  for  tone 
estimation  by  the  complex  subset  selection  method.  The 
examples  are  chosen  from  three  qroups  to  show  the  power  of 
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the  complex  subset  selection  method,  the  windowing 
technique  for  frquency  estimation,  and  the  effect  of  the 
noise. 


In  the  first  example,  the  signal  contains  two  tones 
with  amplitudes  and  frequencies  (A 1= 1. , A2= 1. ) , 

(F1  =  . 2,F2=.4)  ,  respectively.  Twenty  basis  functions  were 
used  where  the  frequency  separation  of  the  basis  functions 
was  .05  (  f=.05  ,  grid  point  separation  ).  Subsets  of  size 
5  (M=5)  were  searched,  and  the  best  one  was  selected  after 

only  examining  136  subsets.  This  is  a  small  fraction  qf  all 
possible  subsets  (20!/(15!5!)  =  15504).  However,  the  two 

frequencies  were  exactly  on  the  qr.id  points.  This  example 
is  illustrated  in  table  4.1  along  with  the  re-ordered  2^ 
errors.  For  the  second  example,  the  amplitudes  and 
frequencies  were  (A  1  =  1 . , A2=1. )  and  ( F 1  — . 2  ,  F2  =  .  585). 
Fifteen  basis  functions  with  frequency  separation  .05  were 
used.  In  this  case,  one  of  the  frequencies  was  off  the  grid 
points,  and  subset  size  8  were  used.  Frequency  FI  was 
selected  exactly  along  with  frequencies  around  F2.  Note 
that  frequency  .6  was  selected  with  a  large  magnitude  which 
is  the  closest  frequency  from  the  set  to  F2.  The  subsets 
examined  were  1716  which  are  less  than  2'  per  cent  cf  the 
total  possible  subsets  (15!/(8!7!)  =  6435).  This  example  is 
illustrated  in  table  4.2  . 

The  third  example  was  picked  to  show  the 
improvement  obtained  by  windowing  tichniquo.  Thus,  the 
amplitudes  and  frequencies  were  chosen  to  be  -the  same  as 


50 


example  2.  In  this  case,  the  selected  frequencies  were  .  2 
and  .  58437  after  examining  only  271  subsets.  It  was 
observed  that  the  windowing  technique  offers  better 
accuracy  while  the  number  of  subsets  examined  are  kept 
small.  This  example  is  illustrated  in  table  4.3  along  with 
three  more  examples  in  table  4.4  to  show  the  efficiency  of 
this  techique.  Note  that  the  results  are  very  accurate, 
even  for  two  tones  with  small  frequency  separation. 

The  next  example  is  an  illustration  of  frequency 
estimation  in  the  presence  of  Gaussian  noise.  Five 
different  signal-to-noise  ratios  (SNR)  were  used,  and  in 
each  case  the  frequency  estimates  were  averaged  over  100 
runs  to  evaluate  the  mean  and  the  standard  deviations.  The 
results  were  tabulated  in  table  4.5  along  with  their 
ccrespondir.g  C-R  bounds.  For  high  SNR'S,  the  standard 
deviations  of  the  estimates  were  consistan*.  with  the  C-R 
bounds,  but  they  differ  for  low  SNR's  due  to  the  occurrc-nce 
of  the  outliers.  Outliers  are  extreme  observations,  and 
they  create  great  difficulty.  When  we  encounter  one,  our 
first  suspicion  is  that  the  observation  resulted  from  a 
mistake  or  other  extraneous  effect,  and  hence  should  be 
discarded.  A  major  reason  for  discarding  it  is  that  under 
the  least  squares  method,  a  fit  curve  is  pulled 
disproportionatly  toward  an  outlying  observation  because  of 
the  squared  diviations  is  minimized.  This  could  cause  a 
misleading  fit  if  indeed  the  outlier  obsevation  resulted 
fi  ud  a  mistake  or  other  extraneous  cause.  On  the  other 
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hand,  outliers  may  con,,cy  significant  information,  as  when 
an  outlier  occurs  because  of  an  interaction  with  another 
independent  variable  omitted  from  the  model. 

Two  frequency  estimation  has  recently  been 
investigated  by  Kumaresan  (15-)  in  a  different  fashion.  In 
his  study,  Kumaresan  searches  over  the  entire  error  surface 
to  locate  the  minimum  error,  and  therefore,  the 
corresponding  two  frequencies.  The  same  example  has  been 
run  by  his  method  and  the  results  are  listed  in  table  4.6 
for  comparison.  Both  methods  give  relatively  identical 
results,  and  consistant  with  C-R  bounds. 

A  small  example  with  10  runs  is  listed  in  table  4.7 
to  show  how  the  output  frequencies  differ  from  the  true 
values  for  each  run.  Two  more  examples  are  illustrated  in 
tables  4.8  and  4.9  to  show  the  estimation  behavior  of  the 
two  frequencies  when  the  frequency  separation  is  very 
small.  The  results  in  table  4.8  are  cor.sistant  wg.th  the  C-R 
bounds  whereas  the  results  in  table  4.9  for  stronger  tone 
are  lower  than  the  C-R  bounds  due  to  the  use  of  information 
about  the  tones  magnitudes.  That  is,  we  have  suplied  more 
information  than  assumed  in  the  C-R  bounds  for  the 
estimation  of  the  frequencies. 

The  C-R  bounds  are  for  unbiased  estimation.  We 
assume  that  our  method,  too,  gives  unbiased  estimation 
since  our  results  are  consistant  with  C-R  bounds,  and  try 
to  show  that  this  assumption  is  correct  by  analyzing  the 
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le  a  ssume,  as  a  working  hypothesis,  that  the 
frequency  error  has  a  normal  distribution  with  zero  mean 
and  standard  deviation  SDT,  Then  if  our  frequency  error 
mean  falls  in  the  interval  -SDT  to  SDT,  we  could  say  that 
the  data  is  consistant  with  our  zero-mean  hypothesis  to 
that  accuracy.  Since  we  have  100  independent  runs,  the 
probability  density  function  of  the  frequency  error  is  a 
chi-sguare  distribution  with  100  degrees  of  freedom  which 
can  be  approximated  accurately  by  a  normal  distibution  (16^ 
having  the  same  mean  and  variance.  Let  us  assume  that  our 

i 

experimental  standard  deviation  (SD)  is  in  error  by  * 20 % 
from  the  true  SDT  (i.e.  SDT  =  (1+,.2)SB).  And  let  us  choose 
the  worst  case,  namely  SDT  -  (1-.2)SD  .  If  the  frequency 
error  mean  is  in  the  interval  -SDT  to  SDT  ,  we  will  say 
that  with  high  probability  ix  is  an  unbiased  estimation.  By 
checking  table  4.5,  we  observe  that  even  for  low  SNR  ,  the 
frequency  error  mean  falls  in  that  interval.  For  instance, 
when  SNR  =  5db; 

FI  -  FI  =  . 906E-3 
SD1  =  . 772E- 2 

SDT  =  ( 1- . 2) SD 1  =  • 6176E-2 
SDT/ (FI  -  *fl)  =  6.82 
-SDT  <  FI  -  ^1  <  SDT 

and  therefore,  our  assumption  is  probably  correct. 
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7ABLL  4.1 

NO.  OF  SAMMLES  =  8 


A 1=0 . 1000D  01  A2= 

■0.10000 

01 

F1-0.2000D  00  F2= 

0.40000 

00 

N=20 

E<  3>=0. 170359 ID- 13 

NO . 

1 

SMALLEST 

ERROR 

E  <  4) =0.2899026D-12 

NO. 

'? 

SMALLEST 

ERROR 

E <  2) =0.57240640-12 

NO. 

3 

SMALLEST 

ERROR 

E <  5) =0.8619378D- 12 

NO. 

4 

SMALLEST 

ERROR 

E (  6 ) =0 . 20045210-11 

NO. 

SMALLEST 

ERROR 

E  <  1 ) =0.3271099D -11 

NO. 

6 

SMALLEST 

ERROR 

E (  7) =0.89958820- 11 

NO. 

7 

SMALLEST 

ERROR 

E  <  8) =0.29528680-10 

NO. 

8 

SMALLEST 

ERROR 

E  <  9) =0.971 80230-10 

NO. 

9 

SMALLEST 

ERROR 

E  <  20  )=().  16025800-09 

NO . 

10 

SMALLEST 

ERROR 

E< 10) =0. 16095810-09 

NO. 

.11 

SMALLEST 

ERROR 

E < 1 1 ) =0 . 1 6555430-09 

NO. 

12 

SMALLEST 

ERROR 

E ( 1 2 ) =0 . 1 6627080-09 

NO. 

13 

SMALLEST 

ERROR 

E< 13) =0.1 6642440-09 

NO. 

14 

SMALLEST 

ERROR 

E < 14 ) =0 . 16646240-09 

NO. 

15 

SMALLEST 

ERROR 

E  <  1 5  )  =0  .  .1 6647270-09 

NO. 

16 

SMALL EST 

ERROR 

E  < 16) =0.16647430-09 

NO  . 

17 

SMALLEST 

ERROR 

£  (18)  =0 . 1 664  7480-  09 

NO . 

18 

SMALLEST 

ERROR 

E< 19)  =0. 16648150-09 

NO, 

19 

SMALLEST 

ERROR 

E<17)=0. 16648690-09 

NO. 

20 

SMALLEST 

ERROR 

N  =  20  M  =  5 


NO.  OF  SUBBETS  EXAMINED  =  1 36 


D< 

1  )  =0.20000 

00 

c< 

1  )  =  (0.1 0000 

0  1  V 

0.76640- 15 

DC 

2)  =0.40000 

00 

c< 

t’J 

II 

*-\ 

c 

:~N 

w 

0.1  . 

0. 16 280-1 4 

) 

D< 

3 >=0.80000 

00 

C( 

3 ) = ( - . 78600 

-  1.  4  y 

-  .  79460  -1  4 

) 

D< 

4) -0.85000 

00 

c:< 

4  )  =  <  0 . 1 .1  660 

13  • 

.101  70-1  3 

) 

D  < 

5) =0.90000 

00 

c:< 

5 >=<0.28960 

-14  y 

0. 10490- 13 

) 

ERROR--  0.90730- 18 


5 

I 

> 

} 

| 

i 

* 


TABLE  4.2 


NO.  OF  SAMPLES  =  10 
A1  =0.:''000D  01 

A2  =o.;  good  01 
FI  =0 ♦  2  good  00 

F2  =0 .  525000  00 

E (  8)-0. 11504850-05  NO.  1  SMALLEST  ERROR 

E<  10)  =0.1 15247911-05  NO.  2  SMALLEST  ERROR 

E (  6) =0.11542740-05  NO.  3  SMALLEST  ERROR 

E (  7) =0.11708150-05  NO.  4  SMALLEST  ERROR 

E <  9 )=0. 11796240-05  NO.  5  SMALLEST  ERROR 

E<13>=0. 12383540-05  NO.  6  SMALLEST  ERROR 
E (  1)=0. 12764290-05  NO.  7  SMALLEST  ERROR 
E< 15) =0. 12934290-05  NO.  8  SMALLEST  ERROR 
E< 14 >=0.12986900-05  NO.  9  SMALLEST  ERROR 
E (  2 >=0.13273340-05  NO.  10  SMALLEST  ERROR 
E <  5 >=0. 13539650-05  NO.  11  SMALLEST  ERROR 
E  <  3>=0. 14358130-05  NO.  12  SMALLEST  ERROR 
E< 11 >=0.15689950-05  NO.  13  SMALLEST  ERROR 
E( 12>=0. 18613920-05  NO.  14  SMALLEST  ERROR 
E(  4  >=0.19928:180-05  NO,  15  SMALLEST  ERROR 
ERROR--  0.183030-05 

D(  1 >=0.20000  00  C  =  <0.999970  00  -.245320-03) 

D(  2  >=0.40000  00  C  =  (-  .442930-02  0 . 777220-02  > 

DC  3 >=0.45000  00  C  =  <-, 366760-  01  -.293640-01) 

DC  4 >=0.50000  00  C  =  (0.105430  00  -.959770-01) 

D<  5 >=0.55000  00  C  =  (0.227950  00  0. 346960  00) 

D<  6 >=0.60000  00  C  =  (0.696760  00  -.314570  00) 

D<  7 >=0.65000  00  C  =  (0.229990-01  0.841350-01) 
D<  8)=0. 70000  00  C,  =  (-.120460-01  0.132500-02) 
NO.  OF  SUBSETS  EXAMINED  =  1716 


TABLE  4.3 


A1=1,0  A2=i  .0 

Fl=.200  F2= .585 


RUN 

/s 

FI 

A. 

F2 

ERROR 

1 

0.2000 

0.60000 

.  68822E-1 

2 

0.2000 

0.57500 

.  31433E-1 

3 

0.2000 

0.58750 

.  19749E-2 

4 

0.2000 

0.58437 

. 12570E-3 

5 

0,2000 

0.58437 

. 12570E-3 

TOTAL  SUBSETS  EXAMINED  =  271 
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TABLE  4.4 


A 1= 1. 0  A2=1 • 0 

No.  of  sanples  =  10 


Case  I 
F1=.  2250 
F2=. 3125 

Case  II 

F 1 =. 2500 
F2  =  , 72  50 

Case  III 
F1=. 1750 
F2=. 5875 

D  1 

D  2 

Error 

0.  2000 

0. 3000 
. 2 1 94E  0 

0.3000 
0.7000 
•  6  7  6  8  E  0 

0.2000 
0.6000 
. 2376 E  0 

WINDOW  1 

D  1 

D2 

Error 

0.  2500 

0. 3000 
.  3485E-1 

0.2500 
0.7500 
. 1827E  0 

0.  1500 
0.6000 
.  203  7E  0 

WINDOW  2 

D 1 

D2 

Error 

0.  2250 

0.  3000 
.2367E-1 

0.2500 
0.7250 
. 2000E-5 

0.1750 
0.5750 
.  4829  E- 1 

WINDOW  3 

D 1 

D2 

Error 

0.  2250 

0.  3125 
. 20  COE- 5 

0. 2500 
0.7250 
.2000E-5 

0.1750 
0.5875 
. 200 0E- 5 

WINDOW  4 

D  1 

D  2 

Error 

0.  2250 

C. 3125 
. 2000E-5 

STOP 

0.1750 

0.  5875 
.2000E-5 

STOP 


\ 


STOP 


CO  VTCTm  £*4 
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TABLE  4.5 


A  1  =  1 .  A  2= 1 . 

F1=. 25  F2=. 375 

NB0N=100  No.  of  samples  =12 


SN  R 

FT 

• 

Ck 

_ 

SD 1 

/> 

F2 

M  -  fcj 

SD2 

20 

. 24987 

.  125E-3 

. 1 36  E- 2 

1 

. 37497" 

. 314E-4 

.  1  48E-2 

10 

.24945 

. 547E- 3 

. 4  24E- 2 

.37478 

.219E-3 

.456E-2 

5 

.24909 

• 906E-3 

.  7  7  2  E  -  2 

.37466 

. 375E-3 

. 628F-2 

0 

. 24850 

.  150E-2 

. 4  84  E-1 

.40651 

.315E-1 

. 126E-0 

-5 

. 32156 

. 715E-1 

.181  E-0 

.‘6109 

.  186E-0 

.  2  37  E-'O 

NF  siqnal  to  noise  ratio 

D  standard  deviation 

i  A  mean  estimate 

i-Fi  frequency  error 


TABLE  4.6 


A  1=1.  A  2=  1 . 

F 1=. 25  F2=. 375 

NFUN=100  NO.  of  samples  =12 


SNR 

!f1  -  FI 
- 1 

SD  1 

1  -A  , 

f2  -  F2j 

SD2 

1 

C-R  BOUNJ 

20 

.  161E-3 

• 128E-2 

. 460E-4 

.  146E-2 

.  133E-2 

10 

. 25QE-3 

. 573E-2 

.230E-3 

.  550  E-  2 

.  420E-2 

5 

. 230E-2 

.470E-1 

. 290E-2 

.050  E-1 

. 748  E-2 

0 

! 

. 339E-2 

.589E-1 

.  338  E-  1 

.  124E  0 

.  133E-1 

1 

-5 

. 186E-2 

•  1 17E  0 

•  1 86E  0 

. 227E  0 

. 236E-1 

TABLE  4.7 

KN=  12  NRUN-  10 

A1=0 •  1000D  01  A2=0.1000Ci  01 

F1=0.2500D  00  F2=0 . 4750D  00 

SIR  s  0.0  SNR  =0 . 1000E  02  NVAR^O . 1000E  00 

FREQ1  FREQ-ERR0R1  FREQ2  FREQ-ERR0R2 

0 . 242187D  00  0.781259D-02  0.478125D  00  0*3124830-02 

0 . 2468758  00  0.312509D-02  0.487500D  00  0.124998D-01 

0 . 259375D  00  0.937491D-02  0.476562D  00  0.156233D-02 

0 . 253125D  00  0.312491D-02  0.47B125D  00  0.312483D-02 

0 . 250000D  00  0.8940700-07  0.470312D  OC  0.468767D-0 

0.248437D  00  0.156259D-02  0.476562D  00  0 . 1 562330-0 

0.253125D  00  0.312491D-02  0*4781250  00  0.3124830-0 

0 . 24B437D  00  0.156259D-02  0.476562U  00  0.156233D-0 

0.2515620  00  0 . 156241D-02  0.4640620  00  0 . 1093770-01 

0 . 250000P  00  0. 894070D~07  0.492187D  00  0.1718730-01 

MEAN1  =  0 . 250312D  00  HEAN2  *»  0.477812D  00 

MERR0Rl=0.312410D-03  MERR(1R2=0 . 2812330-02 

VARIANCE1=0. 1845700-04  VAR  J  ANCE2--0. 55078  ID-04 

STDIV1  =  0 , 429616D-02  ST0IV2  -  0.742146D-02 


r  j  r  j  ro  ro 
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TABLE  4*8 


A1 =1 ♦ 0  A2=l • 0 

Fl=*5  F2  = • 55 

NRUN=100  NO.  of  samples  =12 


SNR 

|pi  -  FI 

■ 

SHI 

1 

|  F2  -  F2 

SD2 

C~R  BOUND 

20 

* 187E-3 

.263E-2 

* 171E-3 

. 262E-2 

. 272E-2 

10 

.281E-3 

*  852E--2 

. 359E-3 

.859E-2 

. 860E-2 

5 

.718E-3 

» 152E-1 

.453E-3 

. 155E-1 

. 153E-1 

0 

* 193E-1 

.8S7E-1 

.594E-1 

. 1 1  BE  0 

*  272E-1 

-5 

.  420E- 1 

. 167E  0 

♦  32 IE- 1 

. 1 79E  0 

.  4  8  4  E  - 1 
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TABLE  4.9a 


11*1.  A  2= 1 0 . 

F1=.25  F2=.  275 

NRUN=100  No.  of  sarples  =12 


SNR 

FI  |F1  -  Flj 

SD1 

F2 

|F2  -  F2[ 

SD2 

20 

.  25078  .78  IE- 3 

.  643  E- 2 

.27515 

.  156E-3 

•  781  E-3 

10 

.25467  . 467E-2 

.224E-1 

. 27402 

. 984E- 3 

|  • 264E-2 

5 

.26035  . 103E-1 

•428E-1 

.27361 

.  139E-2 

.  342E-2 

0 

.27495  .249E-1 

•  7  25e- 1 

.27343 

.  156E-2 

•  391 E-2 

-5 

.29096  .409E-1 

.  8  87  E- 1 

.27309 

! 

.  190F-2  ! 

1 

. 555E-2 

TABLE  4.9b 

A 1  =  1 .  A  2=1 0. 

F 1=. 25  F2=. 275 

NRUN=100  No.  of  samples  =12 


SNR 

C-E  BD 1 

C-P  BD2 

20 

•  7 17E-2 

.  7  17  E-  3 

10 

.  226E- 1 

.  2  26E- 2 

5 

1 

.403E-1 

. 403E-2 

0 

.  717E-1 

.  7 17  E- 2 

5 


127E+0 


127E-1 
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CONCLOSICH 

The  complex  subset  selection  algorithm  was 
described  in  chapter2.  This  algorithm  was  used  in  section 
3.2  to  design  low-pass  PIE  digital  filters  with  a  small 
number  of  coefficients  using  a  weighted  least  squares  error 
criterion.  Digital  lew-pass  filters  with  uniformly-spaced 
time  samples  were  discussed  in  section  3.3,  and  it  was  seen 
in  section  3.4  that,  for  the  same  time  aperture  and  number 
of  coefficients,  the  complex  subset  selection  method 
results  in  a)  smaller  least  sguares  error,  b)  flatter 
magnitude  response  in  the  passband  region,  c)  sharper 
roll-off  in  the  transition  band,  and  d)  higher  peak 
sidelobe  level. 

Tone  estimation  by  the  complex  subset  selection 
method  was  explained  in  section  4.2,  and  was  extended  to 
frequency  estimation  in  section  4.3  .  The  presence  of 

Gaussian  noise  was  discussed  in  section  4.4  .  The  numerical 
results  were  compared  with  the  C-B  bounds  in  section  4.5, 
where  unbiased  estimation  was  assumed,  and  it  was  shown 
that  the  results  were  consistant  with  the  C-B  bounds. 
Finally,  it  should  be  mentioned  that  this  method  can  be 
applied  to  tulti-f reg jency  estimation  with  a  little 
modification,  and  with  seme  more  computation. 


r.  .-•-.vitffii 
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APPENDIX  A 


In  this  appendix,  the  details  of  the  derivation  of 
digital  filler  synthesis  is  presented. 


D(v)  = 


1 

0 

Don't  care 


v  <  p 

si  <  w  <  s2 
,  otherwise 


(A.1) 


W(w)  = 


1  ,  jw  I  ^  p  and  si  ^  |w  |  ^  s2 


0  ,  otherwise 


M 


(W)  =  YL  h  (n)  .EXP  ( JvT  ) 


n-i 


(A. 2) 


(A. 3) 


Where  T  's  are  selected  from  a  set  {ml  with  n=-q, 
T)  '  > 


Min.  Error  = 
{h(D).In} 


«(W) 


D  (w)  -  H  (w)  dw 


(A. 4) 


Using  orthogonality  principle  it  follows 


Jw  (w)  .  (D  (w)  -  H  (w)  )  .  EXP  (-JwT^)  .  dw  =  0 


for  k-1 


H  (A. 5) 
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/' 


H  (w)  .  H  (w)  •  EXF  (~JwT  ).d«  =  W  (w)  .  D  ( w)  .  EXP  (- JvT  )  .  dw 

H  J  X 

for  k=1,...,M  (A. 6) 


A  = 


B  = 


:  B 

(A. 7) 

a!|  a(2 

a2! 

(A. 8) 

a 

•  •  •  3. 

Ml 

mm  _ 

h  (n)  r  •  • 

(A. 9) 

b 

•  .  *  rb  -J 

(A.  10) 

a..  =  {  V  (w)  .EXE-  (Jv  (T-  -T-  )  )  .dw  =  I  EXP  ( Jw.(T;  -T:  ) )  . 

lJ  /p  v  (o':  *<5* 

♦  EXP  (Jw  (!•  -T  •  )  )  .  dw  +  EXP  (Jw  (T-  -T;  )  ) .  dw 
<p  vJ  ^ ^  J 


dw 


=  (  EXP(-Js1  (T.  -T;  )  )  -EXP(-Js2  (T  -T-  ))  )/J(T.  -T*) 
V  J  v  J  L  J 


♦  (  EXP(Jp(T.-T))  -  RXP(-Jp(T.-T,))  ) /J  (T.  -T; ) 

w  VJ  l  l  V. 


♦  (  EXP  (Js2  (Tf  -T;  )  )  -  EXP  (Js  1  (!,'  -T  ;  ) )  )/J(T. -T.*) 

V  yj  L  V)  V  V 

♦  2  [sin(s2(T;-T;))  -  Sin  (si  (T;  -T.  )  ) 

v  0  v  0 

♦  sin  (P  (Tj^-Tj  )  )  3/(T-L-Tj)  (A.  11) 


=  1 


W  (w)  .D  (W)  .EXP  (-JWTO  .  dw  = 


EXP  (- JwT^ )  .  dw 


Therefore, 


-  EXP(JpT^))  3/<"‘7Ty^=  2Sin(pT^)/T^ 

(A.  12) 


H  =  A  B 


(A. 13) 


L 


(  M 

Min.  Error  =  \  « (w)  B(w)  -  )  h  (n)  .EXP(JvT  )  |  .aw 

J  rT*\  * 

h(n),T  M 

=  J  V  (w)  (D  (w)  -  (n)  »EXP( JwT  )  )  .  D(w)  .  dw 


r\z.  I 


n 


(A. 14) 


Min,  Error 


dw 

(w)  .  EXP  ( JwT  )  )  .aw 


i 
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APPENDIX  B 


FOUR  DEFFERENT  PROGRAMS  ARE  LISTED  BELOW 
WHICH  EACH  ONE  IS  FOR  AN  ESPECIAL  PURPOSE. 
PROGRAM  ONE  AND  TWO  ARE  FOR  DIGITAL  FILTER 
DESIGN  WHILE  PROGRAMS  3  AND  A  ARE  FOR  TONE 
AND  FREQUENCY  ESTIMATION. 


PROGRAM -1 


THIS  IS  THE  PROGRAM  LISTING  FOR  DIGITAL  FILTER 
DESIGN  BY  COMPLEX  SUBSET  SELECTION. 


c 

A  (  I  F  J ) 

A  MATRIX 

c 

COL  < I ) 

B  MATRIX 

c 

SOL  ( I  ) 

COEFFICIENT  MATRIX 

c 

Ed) 

El  ERROR  FUNCTIONS 

c 

AE  ( I  ) 

RE-ORDER  ERRORS 

c 

TN<  I ) 

BASIS  FUNCTION  PARAME 

c 

NQ 

Q  IN  EQUATION  (3.2.1) 

c 

M 

NO.  OF  COEFFICIENTS 

c 

P 

PASSBAND  EDGE 

c 

SI 

LEFT  STOPBAND  EDGE 

c 

•  S2 

RIGHT  STOPBAND  EDGE 

c 

T 

GRID  POINT  SEPARATION 

DIMENSION  Y  d  01 ) 

COMPLEX  U f ER , H ( 30 ) * ET , BRSS 
REAL  E <  30 ) , AE  < 30 ) r  TN  <  30 ) , CT (10) 

INTEGER  IF'  <  30  )  » IND  <  30  ) 

INTEGER  IP1 (30) f IP2(30) 

COMPLEX  A <30, 30)  f  At  (30f30)  ,DC0l..(30)  ,  COL  (30) 
COMPLEX  SOL ( 30 > f SOLI (30) 

DATA  A f COL f  Al/900* ( 0 . ,0. ) , 30*0 . 0 f 900* < 0 . , 0 . ) 
DATA  IND » DC0l../30*0 , 30*0 .0/ 

READ ( 5  »  21 0 ) NQ  »  M , NO , P , S 1 , S2 f 7 
READ  <  0 , 220 ) MARK f  NW , NP 
PI =4 ,0*ATAN( 1 .0) 

NQ1~2*NQ  t 1 

IF ( MARK. EQ. 1)60  TO  5 

DO  1  1=1fNQ1 

TN( I )=( I-NQ-1 )*T 


noon 


CONTINUE 
GO  TO  8 

READ(5»230XCT(I>r 1=1 >NU) 

TH=T/2 
DO  7  1=1 fNP 
DO  7  J=1 » NU) 

TN(I+( J-l )*NP)=CT( J)+(I-2)*TH 

CONTINUE 

NQ1=NP*NW 

WRITE(6»150)NG»M 

WRITE<6»240)P»S1 »S2*T 

P=P*2*PI 

S1=S1*2*F‘I 

S2=S2*2*PI 

EVALUATION  OF  THE  ELEMENTS  OF  A  AND  B  MATRIX t 
EQUATIONS  (3.2.14)  AND  (3.2.15) 

DO  32  1  =  1 » NQ1 
DO  20  J=I »NQ1 
T I  .J=TN  ( I )  -TN  <  J ) 

IF(T.l  J.EQ.O.  >00  TO  10 

A(  I  >  J)=2#<SIN(TI  J*S2)~SINU'I  J#S.t  >+SIN<TI  J*P)  )/TI. 
A(J»I)=A(I»J) 

GO  TO  20 

A ( I , J ) =2# <  S2-S1+P )  +  . 0001 
CONTINUE 

IFCTNCI) .EQ.O. )G0  TO  30 
COL( I >=2*SIN<P*TN( I ) )/TN< I ) 

GO  TO  32 
COL  ( I )  --2*P 
CONTINUE 

CALL  AUER IV ( A » COL  >  H  >  ET 1 rNni ) 

ET1 =2#P--ET1 
WRITE ( 6  f 190  >ET1 

EVALUATION  AND  PLOTTING  OF  THE  REOUENCY  RESPONSE 
FOR  THF  FILTER  USING  ALL  THE  BASIS  FUNCTIONS 


N01---N0  M 
L  =  0 

DO  38  I  1 » N 0 1 
V=CMPLX ( 0 . »  0 . > 

DO  37  J=1 f NQ1 

X=2*PI*L*TN< J)/N0 

V=V+H< J)*CEXP(CMPLX(0. »X) ) 

CONTINUE 

L-L+l 

Y  ( I )  =CABS  <  V ) 

Y<  I  >=20*AL0G(  Y(  I  H  .000000.1  ) 


n  n  n  n  n  n  n 
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38  CONTINUE 
P.X--1  ./NO 

CALL  FPL0TT<Y»N01.0. .DX> 2.5 » -100.0) 

EVALUATION  AND  RE-ORDERING  OF  El  ERROR  FUNCTIONS. 
EQUATIONS  <2.2.9)  AND  <2.2.10) 

NS=NQ1-1 
DO  40  IE-1 » NQ1 

CALL  DELT  <  A » COL  » A1 f  DCOL 1 NQ1 » NS  r IE ) 

CALL  AHERIV<A1 » DCOL » SOL » ET1 .NS) 

ET1=2*P-ET1 
E  < IE ) -ET 1 
40  CONTINUE 

DO  60  N=1 »NQ1 
AE  <  K ) -E  <  K ) 

IP  <K  )-K 

DO  50  LP  -1  •  NQ1 
I F  <  AE  <  K ) . LE . E  <  LP ) ) GO  TO  50 
AE <  K  >  ~E <  LP ) 

IP<K)-LP 
50  CONTINUE 

E  <  I P  <  K ) )  =  1 ♦ OE  6 
WRITE (6, 170) IFXK ) *AE<K) »K 
60  CONTINUE 

PEGINING  OF  THE  SEARCH  FOR  BEST  SUBSET 

NR-NQ1-M 
SO  70  1=1 » NR 
I  NIX  I  )  ~1 
70  CONTINUE 

DO  80  J- 1 f  M 
IND< J+NR)=0 
80  CONTINUE 
IE  =  -1 
MP-NR+1 
BR-'AE  <  MP ) 

KP--1 

KPT  =  1 

MI -HU 

ETAN=1.0E  6 

DO  90  IT=1 f MI 

IF < IT. EG. 1 )G0  TO  82 

KP=  <  KP#  <  NR-HT-2  )  )  /  <  T T - 1  ) 

KPT =KPT  +  KP 

82  DO  86  K=1»KP 

IF  < IT ,EQ . 1 )G0  TO  83 
CALL  COMBO < IND.NQ1 rNR.Nf  ) 

83  CALL  INTO  <  A  » COL  *  A 1 .DCOL.Hf TND. IP* IP1 »NR1 jTNj NS  t IE  > 
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CALL  AHERIVIAl fDCOLfSOLfETI fNS) 

ET1=2*P-ET1 

ETA=ET1 

IF IETA.GT .ETAN)GO  TO  86 

ETAN=ETA 

DO  85  IX=1 » M 

SOLI (IX )=SOL  < IX ) 

IP2<IX)=IP1<IX> 

85  CONTINUE 

86  CONTINUE 
IF<ETAN.LE .DR )GO  TO  92 
MP=MP+1 

BR=AE  <  MP  > 

90  CONTINUE 

92  DO  93  L=1*M 

WRITE  <  6  f  1  SOL  » SOL  1<L.)fTN(I  P2  ( L  >  > 

93  CONTINUE 

WRITE  <  6  f  1 90 ) ETAN 
WRITE <  6 . 200 ) KPT 
C 

C  EVALUATION  AND  PLOTTING  OF  THE  FREQUENCY  RESPONSE 

C  FOR  A  BEST  SUBSET  SELECTED  FILTER 

C 

L~0 

DO  94  1  =  1  r NO  1 
V-CMPLX  <  0 . f  0 . ) 

DO  95  J=3  fM 

X-2*F'I*L*TN<IP2<  J>  )/N0 
V~VFS0L1 ( J>*CEXP(CMPLX<0. fX>  ) 

95  CONTINUE 
L=L4-1 

Y< I >-CABS(V> 

Y  <  I  >-20*AL0G<  Y<  I  )  +  .0000001  ) 

94  CONTINUE 

CALL  FPL OT T ( Y  t N01 .  0  . »DX » 2 .5 f - 100 . 0 > 

100  FORMAT (2I2fF10»0) 

110  FORMAT (4F10.0) 

120  FORMAT ( 5X . 'D<'fI2f')='fE10.2) 

130  FORMAT <7 fSXf 'C< ' fI2f ' >-  ('fE12.4f'  f  ' f  E 1? . 4  f  '  >') 

140  FORMAT  (// .  5X  f 'RSS  FOR  FULL  BASIS  ~  ‘  f  E 1 2 . 4  > 

150  FORMAT ( // f 5X f ' N- '.I?.'  M~'fI2) 

160  FORMAT (/fSXf ' E ( ' f 12 . ' ) = ' . El 2 . 4 ) 

170  F0RMAT(/f5Xf  'E<  '  .  I2f  '  )='  fE14. 7f  '  NO.  MP. 

*'  SMALLEST  ERROR') 

180  FORMAT(//f  '  H<  ' fI2f ' ) - (  ' fE1?.4f '  .  '  f f  17 . 4 f ' ) ' f 

*'  T~  ' fE10.4) 

190  FORMAT <//f5Xf 'ERROR  -  ' fL14.7) 

200  FORMAT  < // f  5X  f ' NO .  OF  SUBSETS  EXAMINED  'fI6> 

210  FORMAT (3I3.4F5.0) 

220  FORMAT  <11 .212) 


_ -jWI  AiftaHiri  n*  f'  i’r-  - 


n  n  o  n  n 
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230  FORMAT (5F5.0) 

240  FORMAT ( / . 5X » / PBAND  = ' » E10 . 4 » 3X . ' SBAND1  ='.E10.4.3X, 
♦ ' SBANS2  =' .E10.4.3X. 'T=' ,E10.4) 

STOF 

END 


SUBROUTINE  INTO  IS  FOR  SELECTING  THE  ELEMENTS 
OF  A  ANB  B  MATRICES  FOR  THE  INDICATED  SUBSET 

SUBROUTINE  INTG(A.C0L.A1 . DCQL.H. IND.IP.IP1 »NQ1 ♦  TN.NS. IE ) 
COMPLEX  A (30. 30) »COL < 30) * A1 ( 30.30 ) .  DCOL  <30 ) »H( 30) 

REAL  TN(30) 

INTEGER  IND (  30 ) » IF* ( 30 ) » IPi  (30) 

25  K=0 
K 1  - 1 

DO  26  K2"1»NQ1 
IF ( IND  <  K2 ) . EQ  *  1 ) GO  TO  26 
IPI <K1 )-IP<K2) 

K=K+1 
Kl-Kl+1 

26  CONTINUE 
NS-K 

30  DO  2  1-1  *  NS 
DO  2  J- I » NS 

A1 ( I* J)~A< IPI ( I  >  . IPI ( J) ) 

A1 < J. I )=A1 ( I » J  > 

2  CONTINUE 
DO  3  K-l »NS 
DCOL  <  K ) -COL  < IPI <  K ) ) 

3  CONTINUE 
RETURN 
END 

C 

C 

C  SUBROUTINE  AUER IV  IS  FOR  SOLVING  THE  LINEAR 

C  EQUATIONS  AH  - B»  AND  EVALUATING  THE 

C  APPROPRIATE  ERROR 

C 

SUBROUTINE  AHL'RI  V  <  A  .  COI  .SOL  .LT1  ,K) 

COMPLEX  A  <  30 » 30 )  .COL  (30)  .RL  (30-30)  -SOL  (30)  »E1 .E2.E3.FT 
DO  5  I  -=  1  *  K 
DO  5  J~1 .  K 

5  RL ( I . J ) -CMPLX ( 0 .  .  0 .  ) 

DO  60  1  =  1. K 
1P1  =  IH 
IM1 -1-1 

IF(I .NE .1 > GO  10  30 
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DO  20  J=2 » K 

RL  <  J » 1 ) =CON JG  <  RL  < 1 f  J  > >/RL(l,l> 

GO  TO  60 
DO  40  J=  I » K 
RL(I,J)=A<1» J) 

DO  40  M=1 » IM1 

RL ( I » J )  =RL  <  I » J  >  -RL  <  I » M  )  *RL  <  H » J  > 

IF( I . EQ»K)GO  TO  60 
DO  45  J=IPlfK 

RL  <  J » I ) =CON JG  <  RL  < I »  J ) ) /RL  < I f I ) 

CONTINUE 

DO  100  1=1 fK 

IM1=I-1 

SOL ( I >=COL  <  I ) 

IFd.EQ.DGO  TO  100 
DO  90  L=1 r IM1 

SOL  < I ) =SOL ( I ) -SOL  <  L ) *RL <  I » L  > 

CONTINUE 
DO  140  1=1 fK 
J=K-I+1 

IF  <  J . EQ  »  K ) GO  TO  130 
JF'1  =  JH 

DO  120  L  =  JF‘l  f  K 

SOL  <  .J )  =SOL  <  J )  -SOL  ( L  )  *RL  (  J  f  L  ) 

SOL  <  J ) =SOL  <  J ) /RL  <  J f  J ) 

CONTINUE 

EVALUATION  OF  ERROR »  EQUATION  (3.2.17) 
ET =0 . 

DO  150  I = 1 f K 
ET=ETTCOL(I)*SOL(I) 

CONTINUE 

ET1=CABS<ET> 

RETURN 

END 


SUBROUTINE  COMBO  GENERATES  STRING  OF  ONFS  AND 
ZEROS  FOR  DELETION  OF  BASIS  FUNCTIONS 

SUBROUTINE  COMBO ( JN » NQ1 . M r NF ) 

DIMENSION  JN  <  30 ) 

ONES=M 

NF=0 

NF=1 

CONTINUE 

IF ( JN ( NP ) . EQ . 0 )G0  TO  10 
IF  <NF' . EQ.NQ1  )G0  TO  101 
JN  <  NF* )  =0 


if 
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NF'=NF+1 
ONES-ONES-1 
GO  TO  2 
10  JN<NP >=1 

ONES=ONESH 

IF  <ONES.EQ»M)GO  TO  100 
GO  TO  1 
101  NF=1 
100  RETURN 
END 
C 
C 

C  SUBROUTINE  BELT  IS  FOR  DELETING  ROW  AND 

C  COLUMN  OF  A  AND  B  MATRICES  FOR  EVALUATION 

C  OF  El  ERROR  FUNCTIONS 

C 

SUBROUTINE  BELT  <  A  ,  COL  *  At  *  DCOL.  ,  NG1 , NS , IE ) 
COMPLEX  A <30, 30) ,C0L<30> ,A1 (30,30) , DCOL < 30) 
IEM=IE-1 

IF  < IE . EQ . 1 ) GO  TO  30 
DO  10  1=1,. IEM 
DCOL  < I ) =COL < I ) 

DO  10  J=I , IEM 
A1 < I r J)=A< I » J) 

A1 <J»I >=A< J,I ) 

10  CONTINUE 

DO  20  1=1, IEM 
DO  20  J= IE , NS 
A1 < I f  J ) =A  < I » Jll ) 

A1 < J, I )=A< J+l , I ) 

20  CONTINUE 

IF  <  IE  .  EQ  .  NQ1 )  GO  TO  50 
30  DO  40  I=1E , NS 

DCOL  < I ) =COL  <1  +  1  ) 

DO  40  J=I,NS 

A1  < I »  J ) =A ( I  +  1 , J+l ) 

Al< J,I)=A< J+1,I+1) 

40  CONTINUE 
50  RETURN 
END 


SUBROUTINE  FF’LOTT  IS  FOR  PLOTTING  THE 
FREQUENCY  RESPONSE  OF  THE  FILTER 

SUBROUTINE  FPLOTT < Y » NDIM , XI , DX , YMAX , YMIN ) 

DIMENSION  Y <NDIM ) 

DATA  I BLANK , I POINT , I AST , I MINUS/ '  1PRI 

INTEGER  PLOT (78) 

IF (YMAX-YMIN >250,5,45 


massa 
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5  YMAX=Y(i> 

YMIN=Y  <  1  > 

DO  40  K=1 pNDIM 
IF  <  Y  <  K ) -YMAX ) 20 » 20  p 10 
10  YMAX=Y(K> 

20  IF(Y  <K)-YMIN>30p40p40 

30  YMIN=Y  <  K) 

40  CONTINUE 

IF ( YHAX-YHIN ) 250 p 250 » 45 
45  X=XI-DX 

IF  <YMIN>60p60p50 
50  10=0 

GO  TO  70 

60  IF < YMAX ) 62 p 68, 68 

62  10=76 

GO  TO  70 

68  IQ=IFIX<76,*<-YMIN/<YMAX~YM:[N>  >  +  .5> 

70  K I =10+1 

WR ITE  ( I  PR  I  NT  1 1000  )  YM.IN 
1000  FORMAT  < '  YMIN=  '  p1PE10.3> 

WRITE ( I PRINT p 2000 ) YMAX 
2000  FORMAT  < '  YMAX=  'p1PE10,3/> 

DO  200  K=1 pNDIM 
YK=Y(K> 

IF  (  YK-YMAX ) 81 p  84  p 80 

80  YK=YMAX 
GO  TO  83 

81  IF  <  YK~ YMIN ) 82 p  84  p 84 

82  YK=YMIN 

83  IGUE=IAST 
GO  TO  85 

84  IQUE=IBLANK 

85  CONTINUE 

I P= I F I X  <  76 . *  <  <  YK-YM I N ) / < YMAX - YM I N ) )  +  . 5 ) 

I 6=1 P+1 
X=X+DX 

I <  YK  )  9 0  >  1  50  p 150 
90  DO  100  1=1 p IP 

100  PLOT  < I >=IBLANK 
PLOT ( IB )  =  IAS  T 
J=I B+l 

DO  110  L=J»IQ 
110  PLOT  <  L )  =  IMINUS 

IF(IP-IQ)1 30  p 1 20  p 1 30 
120  PLOKKI  )  =  IAST 

GO  TO  140 

130  PLOT  < KI )  =  1 POINT 

140  WRITE  < I PR I NT  p  3000  )XpY(K)pI  QUE  p  (PLOT  (M)  •  M-- 1  p  K  I  > 

GO  TO  200 

150  DO  160  1=1 pKI 
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160  PLOT ( I )=I BLANK 
PLOT  <  KI )= I POINT 
J=KI+1 

DO  170  I^J»IP 
170  PLOT  < I )=IMINUS 
PLOT ( IB) -I AST 

WRITE ( I PR I NT  f 3000  >  X  *  Y  <  K  >  > IQUE  t < PLOT <  M ) >  M=1 > IB ) 
3000  FORMAT  < '  X~ ' » 1PE10 . 3 r 3X » '  Y= ' , 1PE10 . 3 > 3X , A1 » '  \ 
200  CONTINUE 
RETURN 

250  WRITE ( IPRINT 1 5000 ) 

5000  FORMAT < '  ERROR — YMAX=YMIN ' ) 

RETURN 

END 


» 77A1  ) 
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PROGRAM-2 


THIS  IS  THE  PROGRAM  LISTING  FOR  UNIFORMLY- 
SPACE  PARAMETER  FILTER 


EVALUATION  OF  THE  ELEMENTS  OF  A '  AND  B '  MATRICES 
OF  EQUATION  <3.3,3) 


c 

A( I r J) 

A '  MATRIX 

c 

B  ( I  ) 

B '  MATRIX 

c 

H<  I ) 

COEFFICIENT  MATRIX 

c 

M 

NO.  OF  COEFFICIENTS 

c 

P 

PASSBAND  EDGE 

c 

SI 

LEFT  STOPBAND  EDGE 

c 

S2 

RIGHT  STOPBAND  EDGE 

c 

ND 

NO.  OF  DELTAS  FOR  S 

c 

NC 

NO.  OF  C'S  FOR  SEAR 

c 

DEL 

INITIAL  DELTA 

c 

DELD 

DELTA  INCREMENT 

c 

Cl 

INITIAL  C 

c 

DELC 

C  INCREMENT 

DIMENSION  DN<10) >Y<101 ) >E<100) 

COMPLEX  H<10>»A(10»10>»B<10> »HB<10> »ET»V 

READ  <  5  > 1 00 ) NO » M » ND  r  NC 

READ  <  5  r 1 1 0 ) P . SI , S2 » DEL » DELD 

READ ( 5  f 1 10) Cl *DELC 

PI=4*ATAN< 1 .  ) 

D1=DEL 

P=2*PI*F‘ 

S1=2*PI*S1 
S2=2*PI*S2 
SER=1.0E  6 
Kl=4 
CO=CI 

DO  50  IK~1 t NC 
DO  40  IL=1 » ND 
DO  25  1  =  1  >M 
DO  15  J=I »M 
IF<T.EG.J> GO  TO  10 
X=D1*< I-J) 

A  < I »  J  >  =2*  <  SIN<  S2*X  > -SIN ( S1#X ) +SIN<  P*X ) )/X 
A< J, I )=A< I , J) 


GO  TO  15 

10  Ad  »  J>  =2* ( S2-S1+P  >  +  » 0001 
15  CONTINUE 

X=d-3)*D1+C0 
IF<X.EQ.O. ) GO  TO  20 
Bd)=2*SIN<P*X)/X 
GO  TO  25 
20  Bd  >=2#P 
25  CONTINUE 

CALL  AHERIV< A *B,H»E1,ET,M> 

SELECTION  OF  THE  SMALLEST  ERROR  AND 
CORRESPONDING  COEFFICIENTS 

E< IL )=E1 
ET1=2*P-E1 

WRITE  <  6 » 125>CQ*  D1 »ET1 

IF ( SER . LT . ET1 >  GO  TO  35 

SER=ET 1 

CB=CO 

BD=D1 

DO  30  L=1,M 
30  HB  <  L ) -H  <  L  > 

35  D1=D1+DELD 
40  CONTINUE 
CO=CO+DEL.C 
D1=DEL 

50  CONTINUE 

WRITE ( 6 » 1 30 ) CD , BD » SER 
DO  60  1=1, M 
60  WF.’  I TE  <  6 , 1 4  0 ) HB  d  ) 

DO  62  1  =  1  *M 
62  DN d  )  =  < 1-3 ) #BD+CB 

WR I TE  <  6 » 1 50 ) <  DN  < I > ,  I  =  1  ,  M ) 

EVALUATION  AND  PLOTTING  OF  FREQUENCY  RESPONSE 

N01=N0H 

L=0 

DO  70  I  =  1 » NO 1 
V=CMPLX(0. »0, ) 

DO  65  J=1 , M 

X=2*PI*L*DN( J>/NO 

V=V+HB<  J>*CEXP<CMF'LX<0.  ,X>  ) 

65  CONTINUE 
L=L  +  1 

Y( I )=CABS< V) 

Y(I) =20*AL0G <  Y  d  )  +  . 0000001 > 

70  CONTINUE 
DX=1 ./NO 
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CALL  FF'LOTT  <  Y » NO!  » 0 . » DX , 4  » , - 100 . 0  > 

CALL  FPLOTT  <  E » ND ,  DEL , DELC ,  1 » 8 , 0 . 0  > 

100  FORMAT  (13,312) 

110  FORMAT (5F5.0) 

120  FORMAT < / , 5X  *  "  DELTA= '  , E 10 . 4  *  '  ERROR=  '»3E12»5> 

125  FORMAT  <  3X ,  'C='  , E9 . 3 , 3X , '  DELTA=  '  » EC! » 2 , 3X ,  ' ERROR-  '  ,  E12 » 6 ) 
130  FORMAT </»5X» 'BEST  C='»E10,4»'  BEST  DELTA= ' , E10 . A , 

*'  £RROR= ' » 2E 1 2 . 5 ) 

140  FORMAT  < / , 5X ,  '  HB  =',2E12.5> 

150  FORMAT  </»5X»5<E10.4,2X> ) 

STOP 

END 


SUBROUTINE  AHERIU  IS  FOR  SOLVING  THE  LINEAR 
EQUATIONS  A ' H -B '  ,  AND  EVALUATING  THE 
CORRESPONDING  ERROR  FUNCTION 

SUBROUTINE  AHERIV ( A , COL , SOL , El , ET » K ) 

COMPLEX  A  < 10 , 1 0  > , COL  < 10  > , RL  < 1 0 , 1 0  > ,  SOL  < 10) ,E1 ,E2,E3, ET 
DO  5  1  =  1, K 
DO  5  J= 1 ,  K 

5  RL ( I » J  > =CMPLX  <  0 .  ,  0  *  > 

DO  60  I  =  1 »  K 
IP1  =  I +1 
IM1=I~1 

IF ( I . NE  *  1 ) GO  TO  30 
DO  10  J=1  » K 
10  RL  <  1 ,  J  >  =A  <  1 ,  J  > 

DO  20  J=2»K 

20  RL ( J , 1 ) =CON JG ( RL <  1 ,  J ) ) /RL  (1,1) 

GO  TO  60 
30  DO  40  J=I,K 

RL  < I , J  > =A <  I ,  J  ) 

DO  40  M=1 , IM1 

40  RL  < I , J )  =RL  <  I  »  J  )  -RL  <  I  »  M  >  *RL  <  M  •  ,J ) 

IF ( I . EG . K )G0  TO  60 
DO  45  J=IP1 *  K 

45  RL  <  J  ,  I ) =CON JG  <  RL  < I  , J ) ) /RL (1,1) 

60  CONTINUE 
80  DO  100  1=1, K 
IM1=I-1 
SOL  ( I )  =COl.  ( I ) 

IF  < I . EQ . 1 )G0  TO  100 
DO  90  L=1 , IM1 

90  SOL < I ) =SOL < I ) -SOL ( L ) *RL ( I , !  ) 

100  CONTINUE 

DO  140  1=1, K 
J=K-I+1 

IF  <  J , EQ  >K ) GO  TO  130 
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IF<J.EQ.K)GO  TO  130 
JP1=J+1 

DO  120  L=JF‘l  »K 

120  SOL  <  J ) =S0L ( J ) -SOL  <  L ) #RL  <  J » L  > 

130  SOL  <  J )  =SOL  ( J )  /RL  <  J » J ) 

140  CONTINUE 

EVALUATION  OF  ERROR »  EQUATION  <3.3.4) 
ET=0. 

DO  150  1=1 »K 
ET=ET+COL< I )#SOL< I ) 

150  CONTINUE 

E1=CABS<ET> 

RETURN 

END 
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c 

c 

c 

c 


PROGRAM-3 


c 

THIS  IS 

THE  PROGRAM  LISTING  FOR 

c 

c 

TONE  ESTIMATION. 

c 

c 

A1  >A2 

TONES  AMPLITUDES 

c 

B1  .  B2 

-ONES  FREQUENCIES 

c 

KN 

,;0.  OF  SAMPLES 

c 

A  ( I » J  > 

A  MATRIX 

c 

COL ( I ) 

B  MATRIX 

c 

SOL  <  T ) 

COEFFICIENT  MATRIX 

c 

E(I) 

El'S  ERROR  FUNCTIONS 

c 

SIR 

SIGNAL  TO  INTERFERENCI 

c 

SNR 

SIGNAL  TO  NOISE  RATIO 

COMMON  PlfAl f A2.B1 1 B2.D1  . D2r D3.D4.N1 r  N2  > KN 
REAL*3  D ( 50 ) » E  < 50 ) » AE ( 50 ) .  A 1 » A2 » B 1 . B2 » P I ,  D 1 » D2 
REAL *8  03  *  D4 , PI » P2 , EC 3 >  P3  * BR . ETA  *  E T AM 
INTEGER  IP (50) , IND ( 50 ) » INBB < 50 ) * INP(SO) > JP(50) 


INTEGER  JJP  (50),  l.JP  <  50  >  .  R 

C0MF'LEX*1  A  A  ( 50 .  50 )  *  COL  <  50  )  ..  SOL  ( 50  )  >  SOI.  1  <  50  )  .  ET 
COMPLEX* 16  S(50>50) ,DC0L(50) » BRSSr Y ( 100) .UN.E3 
DATA  A . COL . I ND/2500* ( 0 . ODO . 0 . ODO ) , 50* ( 0 . DO . 0 , DO ) . 50*0/ 
PI -4  *ODO*DATAN( 1 .ODO) 

READ (5.110) SNR . S I R . P 1 . B2 

IND 1-0 


IND3-0 

A1 =1 .DO 

A2 =A1*( 10** < -SIR/20) ) 

READ (5.1 05 ) KN . AK 
IX=3A57137 

SD=A1*( 10** ( -SNR/20) ) 
VAR--5Sn**2 
WRITE (6,1 13)KN 
WRITE ( A . 1 1 5 )  A1  . A2 
WRITE  (  A  »  1 1 7  )  PI  fh'2 
10  READ  (  5 . 100  >  L.F  I  »  N . MS . ML . D1  .  D2 
IECLFI .EG. 1 )G0  TO  2 


GENERATION  OF  BASIS  FUNCTION  PAR AMR  If RS 


non  o o  o  n 
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n<i >=di 

DO  3  1=2, N1 
3  D< I >=D< I-i )+D2 

D(N1+1 )=D3 
DO  4  J=2,N2 

A  D(N1+ J)=D(N1+J— 1 )+D4 

N=N1+N2 

5  WRITE(6,210)N 

I F ( I ND 1 » EG . 1 ) GO  TO  11 

GENARATION  OF  GAUSSIAN  NOISE  IN  DATA, 

EQUATIONS  (4.2.2)  AND  (4.4.1) 

DO  7  1=1, KN 

CALL  GAUSS (IX, SD , 0 . , V  ) 

CALL  GAUSS  (IX,  SD  ,  0  .  ,  V.1  ) 

VN  =CMF'LX  ( V » VI  ) 

P1=2.D0*PI*B1 

P2=2.D0*PI*B2 

Y(I  )=A1*CDEXP(DCMPLX(0.D0,P1*(  1-1  )  )  ) 

Y ( I )=Y ( I  >+A2#CDEXP( DCMPLX  ( 0 .  DO ,  F'2#  ( I  - 1 )  )  > 

7  CONTINUE 

11  IF ( IND3 . EQ . 0 >GQ  TO  8 
P1=2.*PI*D(LJP(1 ) ) 

P2=2.*PI*D(LJP<2) ) 

DO  9  1=1 ,KN 

Y(  I  )=Y(I  )-SOLl  (1  >*CDEXP(DCMPLX(O.DO,Pl*(  X-l  )  )  ) 
Y ( I ) =  v ( I ) -SOL 1(2) *CDEXP ( DCMPLX ( 0 . DO , P2* < I  3 ) ) ) 
9  CONTINUE 
0  NS~N 

EVALUATION  AND  RE-ORDERING  OF  El  ERROR 
FUNCTIONS,  EQUATIONS  (2.2.9)  AND  (2.2.10) 

IE  =  0 

CALL  I NTG ( A , COL , Y , I ND , IP , N . NS , I  F I , It  ) 

CALL  AHER I V ( A , COL , SOL , Y , E T , NS , KN  ) 

ETA=CDABS ( ET ) 

NS=N-1 

DO  40  IE=1,N 

CALL  BELT ( A , COL , S , DCOL ,N  » NS , IE ) 

CALL  AUER I V ( S , DCOL , SOL , Y , ET , NS , KN ) 

E ( IE ) =CDABS ( ET ) 

40  CONTINUE 

DO  60  K=1,N 
AE(K)=E(K) 

IP(K)=K 
DO  50  LP=1,N 

IF(AE(K) .LE.E(LP) )G0  TO  50 
AE(K)=E(LP> 


o  n  ft 


IP  <  K ) -LP 
50  CONTINUE 

E  <  IP  <  K ) )  =  1  *  OE  6 
60  CONTINUE 

DO  2000  IXX-1 »N 

WRITE  <  6  > 170>IP<IXX) » AE  <IXX)»IXX 
2000  CONTINUE 

BEGINING  OF  SEARCH  FOR  THE  BEST  SUBSET 

BO  99  M=MS  > ML 
WRITE (6*1 50 )NtM 
NR=N~M 

no  70  1=1 » NR 
iNri<D=i 
INDB  < I ) =1 
70  CONTINUE 

BO  80  J=i,M 
INB( J+NR  >=0 
INDB<  .JiNR)=0 
80  CONTINUE 
IE--1 
MP=NR+1 
BR=AE  <  MP  > 

KP=  1 

KPT  =  1 

M I  =  M  M 

ETAN=1 . OE  6 

BO  90  IT -1 »  MI 

IF  < IT . EG . 1 )G0  TO  82 

KP=  <  KP*  <  NR+IT-2  > )/< IT-1 ) 

KPT =KPT+KP 

82  BO  86  K  =  1  t  KF' 

IF < IT . EG *  1 ) GO  TO  83 
CALL  COMBO  < INB *  N  *  NR  >  NF ) 

83  CALL  INtG  <  A » COL  f  Y  t  INB » IP r  N » NS  » LF I  »  Il-Z ) 
CALL  AHER I  V  ( A » COL  »  SOL  » Y » ET  » NS » KN  > 
ETA=CHABS(ET> 

IF  <  ETA . GT . ETAN  > GO  TO  86 

ETAN=ETA 

ETAN -ETA 

BO  84  1-1 »N 

INBB  ( I )  =  INB  < I  ) 

84  CONTINUE 
DO  85 

SOLI < J)=SOl  < J) 

85  CONTINUE 

86  CONTINUE 

87  IF  <  ETAN  .  L.E  .  BR )  GO  TO  92 

88  MP-MP+1 
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BR=AE<MP) 

90  CONTINUE 

92  J=0 

DO  93  L=1 t N 

IF<INDB<L).EG.1)G0  TO  93 
J=J+1 

JF‘<  J )  =  IP  ( L  ) 

93  CONTINUE 

DO  95  l..F'=l  i  J 
L  JP  ( LP )  =  JF‘  ( LF* ) 

JJP  <  LF' )  =LP 
DO  94  KP=1*J 

IF<L.JP(LP)  .LE.JP(KP)  )G0  TO  94 
L JP  <  LP ) =  JP  <  KP ) 

JJP<LP)=KP 

94  CONTINUE 

JP(  J.JP<LP>  >  =  1F  6 

95  CONTINUE 
WRITE (6, 200) KPT 
DO  96  L.  - 1*  J 

WRI TE  <  6  »  1 80 )  L  » D  <  L  JP  ( L  .  >  )  .•  L  >  SOL .1  (  i  ) 

96  CONTINUE 

WRITE  <  6  »  .1  90 )ETAN 

99  CONTINUE 

READ  < 5  j  240 ) IND2 » IND1 - IND3 
I E  < I ND3 • EO . 1 ) GO  TO  10 
IF  < IND2 . EQ » 1 )G0  TO  101 
GO  TO  10 
101  CONTINUE 

1 00  FORMAT ( 1 1  *  31 2  *  2F 1 0 . 0 > 

105  FORMAT < 1 3 » F 1 0 . 0 ) 

110  FORMAT <4F 10,0) 

113  FORMAT ( / » 5X r 7  NO .  OF  SAMMIES  ^'»I4> 

115  FORMAT  < / f 5X  * 7A1  =  7 »E10.4»5X»  7A2=7 »  E 10.4) 
1 1 7  FORMAT  < / r  5X  » 7 FI = ' »  E 1 0 . 4  » 5X  » 7 F2 - 7 f  E 1 0 . 4 ) 
120  FORMAT  <  5X  t  7  D  < ' r 12 r  7 )  ~  7 »  E 1 2 . 4 ) 

125  FORMAT  < / >  5X > 'SIR-' *E10.4>5X» 7 SNR- ' >E10. 
*  7  VAR~ 7  >  E 1 0 . 4  ) 

130  FORMAT  </»5X»  'C< 7  »  1 2  f ' )  -  <7, El  2.4* 7  *  7 

140  FORMAT <//*5X> 7 RSS  FOR  FULL  BASIS  --7»E12 
150  FORMAT  < // >  5X  t  7 N  =  7  » 1 2  * 7  M=7>12) 

1 60  FORMAT  C  /  ,  5X  *  7  E  (  7  ,  1 2 »  7  )  ■== 7  *  F 1  3 . 7  ) 

170  FORMAT  < / 1 5X  >  7E( 7  *I2»  7 )"' *  E 1 3 . 7 »  7  NO. 

*7  SMALLEST  ERROR7) 

180  FORMAT  < / » 5X  »  D<  7  »I2f 7 >  =  ' »E10.4»5V» 7C<  '  ■ 
#E . 3  0 . 4  »  7  ,  7  »E10.4r  7  )') 

190  FORMAT <//»5X»  'ERROR--7  .E12. 4) 

200  FORMAT < // r 5X r  7 NO .  OF  SUBSETS  EXAMINED  - 
210  FORMAT ( // »5X  r  7  N= 7  » 12 ) 

220  F0RMAT<2I2f2F10.0) 


» E 1 0 . 4  » 5  X  > 


*  7  >  El  2 . 4  >  7 

» E 1 2 . 4  ) 


7  »  I  7  » 
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230  FORMAT <//»5X» 'SUM  OF  COEFFICIENTS  SQUARED  ='»E12.4) 
240  FORMAT (311) 

STOP 

END 

C 

C 

C  SUBROUTINE  INTO  IS  FOR  EVALUATION  OF  ELEMENTS  OF 

C  A  AND  B  MATRICES  FOR  EACH  SUBSET  TO  BE  EXAMINED . 

C  EQUATIONS  <4*2.14)  THRU  (4.2.16) 

C 

SUBROUTINE  INTO ( A  »CGL  » Y » IND. IP. N . NS . LFI » IE) 

COMMON  PI  » A1 » A2 » Bi . B2 » D1 . D2 »  D3  » D4  » Nl » N2 » KN 
COMPLEX* 16  A  <  50 » 50 ) . COL  < 50 ) » Y  < 50 ) 

REAL *8  D ( 50 ) »AD<50> . X » Y1 . Y2 . Y3 . A1 . A2 , B1 .B2.PI 
REAL*8  D1 .D2.D3.D4.AK 
INTEGER  IND ( 50 ) » IP< 50 >  » JP ( 50 ) 

IF  <  LFI . EG . 1 )G0  TO  2 
D ( 1 ) =D 1 
DO  1  I..-2  .N 

1  D < L ) ~D < L-l ) +D2 
GO  TO  S 

2  D  <  1  >  •  D 1 

DO  3  L~2.N1 

3  D(L  >=IHL.-1  )+D2 
D<  Nl  f I ) ~ D 3 

DO  4  LL. 2  >  N2 

4  D  <  Nlf  LL.  >  =•■  (Nl+LL  -1  )  +D4 

5  IF ( IE. EG. 0) GO  TO  30 
IF (IE. EG. -1 )G0  TO  25 
NS-N-1 

DO  20  L  =  1 » N 
IF  <  L , NE , IE )G0  TO  20 
DO  15  M-L  »NS 
15  D (M )=D< M+l ) 

20  CONTINUE 
GO  TO  30 
25  K=0 

Nl  =1 

DO  26  K2 -1 .N 

IF ( IND ( K2 > . EQ . 1 )G0  TO  26 

AD(K1 )=D( IP(K2) ) 

K=K+1 
K1*KH  1 
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JP<K3)=K4 

27  CONTINUE 

AD< JR 'R3> )=1 .OE  6 

28  CONTINUE 
NS=K 

30  DO  60  1=1, NS 
DO  60  J=I »NS 
X=<D< J)-D(I) >*2.0D0*FI 
IF<I-J)40»  50 » 40 
40  CONTINUE 

All, J>=<1 . ODO-CDEXP < DCMPLX  < 0 . DO ,  X*KN ) >  > 

A  < I » J  >  = A  <I»J)/<1. DO-CDEXP <  DCMPLX ( 0 . DO » X ) ) > 

A ( J, I ) =DCON JG  <  A  <  I » J  >  > 

GO  TO  60 

50  A  < I »  J ) =1  . ODO*KN 
60  CONTINUE 

DO  80  K=1»NS 

COL  ( K  >  =  ( 0  .  DO » 0 .  DO  > 

Y3=-2,D0*PI*D<K> 

DO  80  J=1 >KN 

COL ( K ) =  COL ( K ) + Y ( J ) *CDEXP <  DCMPLX <  0 . DO  *  Y3* ( J - 1 >  )  > 
80  CONTINUE 
RETURN 
END 
C 
C 

SUBROUTINE  AHERIU  IS  FOR  SOLOING  THE  LINEAR 
EQUATIONS  AC=B  AND  CORRESPONDING  ERROR'  I  UNCTION 

SUBROUT I NE  AHER I U  <  A  »  COL , SOL , Y , ET , K ,  KN ) 

COMMON  PI  » A.l  >  A2 »  B.t  ,B2 ,  D1  ,  D2»  D.3  »  D4  >  N1  »N2 
REAL*S  A1  »  A2  , B1 ,B2,PI ,C,VN<50> rCl »C2 
COMPLEX*  16  A  ( 50 » 50  )  » COL  (50)  ,RI  (50,50)  , SOL (50) 
COMPLEX*  1 6  XXI  ,Y< 100)  ,E.t  ,E2,E3,ET,XX 
DO  5  I  -- 1  ,  K 
DO  5  ,J~  1  *  K 

5  RL  < I >  J )= DCMPLX <  O . ODO *  0 . ODO ) 

DO  60  1  =  1,  N 
IPl=Iit 
I M 1  =  I  ~  1 

IF ( T . NE . 1 ) GO  TO  30 
DO  10  J= 1 »  K 
10  RL  <  1  ,  J  >=A <  1  ,  J  ) 

DO  20  J=2  *  K 

20  RL  <  J  ,  1  )  =DCON JG  <  RL  ( 1 »  J)  )  /RL  (  !  ,  .1  ) 

GO  TO  60 
30  DO  40  J= I  » K 

RL  < I »  J ) =A ( I  ,  J) 

DO  40  M=1 »  IM1 

40  RL  <  I  ,  J  )  =RL  <  I  »  J  )  -RL  (  I  »  M  )  *RL  (Mr.l) 
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IF < I . EQ  * K ) GO  TO  60 
00  45  J=IP1,K 

45  RL ( J , I ) =DCON JG ( RL ( I »  J ) )/RL(IrI) 

60  CONTINUE 
60  DO  100  1  =  1, K 
I M 1  =  I  —  1 
SOL  < I >=COL ( I ) 

IF( I »EQ. 1 )G0  TO  100 
DO  90  L=1 » IM1 

90  SOL ( I > =SOL  < I ) -SOL  <  L ) *RL  < I , L ) 

100  CONTINUE 

DO  140  1=1, K 
J=K— 1+1 

IF(J.EG.K) GO  TO  130 
JF'1=J+1 

DO  120  L--.JP1  ,K 

1 20  SOL  <  J )  =SOl.  <  J  )  -SOL  <  L  )  *RL  <  J  ,  L  > 

130  SOL  (  J ) -SOL ( J >/RL  <  J, J) 

140  CONTINUE 
C 

C  EVALUATION  OF  ERROR,  EQUATION  (4.2.19) 

C 

E2=0 . ODO 
E 1 =0 . ODO 
DO  15  11=1, K 

E2=E2iS0L  (ID  *DCON  JG  <  COL  (ID) 

15  CONTINUE 

DO  16  1 3= 1 , KN 
E 1 =E 1 1 C DABS ( Y ( 1 3 ) >  **2 

16  CONTINUE 
ET=(E1  ~E2)./KN 
RETURN 

END 

c 

r: 

SUBROUTINE  COMBO  IS  FOR  GENERATING  STRING 
OF  ONES  ANI  ZEROS  FOR  DELETION  OF.  PROPER 
0  BASTS  TUNC  JONS 

C 

SUBROIJT I  M.  COMBO  (  JN  ,  N  ,  M ,  NF  ) 

DIMENSION  N .40) 

IONES=M 

NF---0 

1  NP  =  1 

2  CONTINUE 

IF(  JN’(NF')  .  t  .0)00  7  i  1<> 

IF  'NE'.EO.  N>  0  TO  101 
JN ( NP ) =0 
NF--NPU 
1 0NE5  =  I  ON!  S  1 
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GO  TO  2 
10  JN<NP) =1 

I0NES~I0NES+1 
IF  < I ONES ► EG  *  h) GO  TO  100 
GO  TO  1 
101  NF=1 
100  RETURN 
END 
C 
C 

C  SUBROUTINE  GAUSS  AND  RANDU  ARE  FOR 

C  GENERATING  GAUSSIAN  NOISE 

C 

SUBROUTINE  GA'JSSdX.S.AM.V) 

A~0 * 0 

DO  50  I  —  1 y 1 2 

CALL  RANDU (IX.IY.Y) 

I X  ~  I Y 
50  A-A1Y 

V=  <  A-6 ► 0 ) #S+AM 

RETURN 

END 

SUBROUTINE  RANDU < I X  > I Y  >  YFL  > 

IY=IX*65539 

IF<IY)5»6>6 

5  I Y= I Y+ 2 14748364 7+1 

6  YFL --I Y 

YFL-YFt *.4656613L~? 

RETURN 

END 

C 

C 

C  SUBROUTINE  BELT  IS  FOR  DELETING  F'ROFFR  ROW 

C  AND  COLUMN  OF  A  AND  B  MATRICES  FOR 

C  EVALUATION  OF  El  ERROR  FUNCTIONS 

C, 

SUBROUT  I NE  BELT  <  A  >  CUL  > S  »  DCOL. >  N , NS  > I E  > 

COMPL.  I  X*  I  6  A  (  50  »  50  )  >  COL  ( 50  )  >  S  (  50  >50)  *  DCOL  (  50  ) 
TEM-IF  -1 

IF  <  IE  .  E'O  *  1 )  GO  TO  30 

nn  10  it.  tem 

DCOL < I ) -COL < I ) 

DO  10  J~I . IEM 
S  < I >  J ) ~A  < I »  J ) 

S  ( .1.  I )  A  ( ,  > »  I  ) 

1 0  CON  r I nul 

no  20  1-1 f JEM 
DO  20  J= IE. NS 
S<  I »  J)=A<  I  ,.J+1  ) 

S< J. I )=A<  H 1 .  I  ) 
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20  CONTINUE 

IF<IE.EQ.N)GO  TO  50 
30  DO  40  I=IE»NS 

nCOL(I)=COL(I+l) 
no  40  J=I,NS 
S  < I »  J )  =  A  < 1  +  1 »  J+l ) 

S< J>I)=A< J+1»I+1) 

40  CONTINUE 
50  RETURN 
ENEi 
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C  PROGRAM-4 

C 
C 
C 

C  THIS  IS  THE  PROGRAM  LISTING  FOR  TWO 

C  FREQUENCY  ESTIMATION  BY  COMPLEX  SUBSET 

C  SELECTION  AND  WINDOWING  TECHNIQUE. 

C 
C 

C  A1.A2 

C  B1 *  B2 

C  KN 

C  A  <  I » J  ) 

C  COL < I ) 

C  SOL  < I ) 

C  Ed) 

C  SIR 

C  SNR 

C  NRUN 

C 

COMMON  P I » A 1 >  A2 . B 1 . B2 . D 1 . D2 »  D3 .  D4  » N1 . N2 » KN 
REAL *0  D <  40 ) .E <  40 ) . AE  <  40 > . AJ . A2. B1 . B2 . PI 
REALMS  Dt .D2.D3.D4.BR.ETA.ETAN.pl .P2 
REALMS  EDI 1 <100) »ED22<100> » DMER1 >  DMER2  >  STD1 . STD2 
REAL  *B  EDI <100). ED2 (100). DMEAN 1 . DMCAK'2 . OAR t . 0AR2 
INTEGER  1 1 '  <  4  0  )  »  IND<40)  ,  INDB<40)  .  INP  <  40 )  .  JP  <  40  ) 
INTEGER  L  JP  <  4  0 ) . R . J JP (40) 

C0MPLEX*16  AMO. 40)  .C0|„<40>  .S0L<40)  .SOI  1  <40)  .ET 
COMPLEX* 16  UN. S< 40. 40) .DCOL  <40) , BRSS . Y < 1 00 ) 

DATA  A  >  COL/ 1 600*  <  0 . ODO . 0 . 0 DO ) . 404 <  0 . DO . 0 . DO ) / 
DATA  1 ND/40*0/ 

P I “4 . ODO*DAT  AN  < 1 .ODO) 

READ <5.110) NRUN .SIR. SNR . B 1 . B2 
READ  <  5 . 1 05  )  KN » N1  .  MS  >  Ml. 

A 1=1 .DO 

A2=A1*  < 10** < -SIR/20 ) ) 

WR I TE  <6.230) KN . NRUN 
WRITE  <6.35<))A1  »A2 
WRITE <  6 . 360 >B1 .B2 
KPT--0 
1X^3657137 

SD=A1*  < 10**  <  —SNR/20 ) ) 

VAR=SD**2 

WR I TE  <  6 . 280 ) S I R . SNR . UAR 
Pl=2.D0*fI*Bl 
P2=2.D0*PI*B2 
C 

C  BEG IN. TNG  OF  EACH  RUN 

C 

DO  500  I J- 1 .NRUN 
C 


TONES  AMPLITUDES 
TONES  FREQUENCIES 
NO.  OF  SAMPLES 
A  MATRIX 
B  MATRIX 

COEFFICIENT  MATRIX 
El  ERROR  FUNCTIONS 
SIGNAL  TO  INTERFERENCE  RATIO 
SIGNAL  TO  NOISE  RATIO 
NO.  OF  RUNS 


oono  nno  non 
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GENERATION  OF  GAUSSIAN  NOISE  IN  DATA* 
EQUATIONS  <4.2.2)  AND  <4.4.1) 

DO  1  1  =  1  *  KN 

CALL  GAUSS (IX*  SD  *  0 . » V ) 

CALL  GAUSS  < IX  *  SD  *  0 . *  U1 > 

UN  CMPLX < V* U1 ) 

Y  ( I )  =A1 *CDEXP<  DCMPLX  <  0 .  DO  *  F*1  *  <  I  ~1 )  )  ) 

Y  < I ) =Y ( I ) +A2*CDEXP  <  DCHPLX  <  0 . DO  * P2*  <  I  ■- 1  )  >  )  +  UN 

1  CONTINUE 

GENERATION  OF  BASIS  FUNCTION  PARAMETERS 

N  =  10 

N0=0 

LF  1  =  0 

D 1  =  .  1 

D2=D1 

D  ( 1  >  --0 .  DO 

DO  2  .1=2*10 

2  D< I >=D< 1-1 )+Dl 
NS=N 

GO  TO  7 

3  LFI=1 
D< 1 >=D3 
D  <  A ) =D4 

DO  4  J=2  *  5 
D<  J)=D<  J-D+Dl 

4  D<  J+5)=D<  J+4M-D1 
N'-l  0 

NS-N 
7  IE=0 

CALL  INTO  <  A  *  COL  *  Y  * IND , IP  *  N  *  NS  *  LFI *  IE ) 

CALL  AUER I U  <  A  *  COL  *  SOL  *  Y » ET  *  NS  * KN ) 

ETA=CDABS  <  ET  > 

EUALUATION  AND  RE-ORDERING  OF  El  ERROR 
FUNCTIONS*  EQUATIONS  <2.2.9)  AND  <2.2.10) 

NS=N-1 

DO  40  IE=1*N 

CALL  BELT  <  A  *  COL  *  S  *  DCOL  *  N  *  NS  *  I E ) 

CALL  AHER I U  <  S  *  DCOL  *  SOL  *  Y  *  ET  *  NS  *  KN ) 
E<IE)=CDAPS<ET> 

40  CONTINUE 

DO  60  K  - 1  *  N 
AE  <  K ) =E  <  K ) 

IP  < K  ) ~-K 
DO  50  LP--1  *  N 

IF  <  AE  <  K ) , LE . E <  LP ) ) GO  TO  50 
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AE<K)=E<LF) 

IF*  (K )  =LP 
50  CONTINUE 

E  < IP < K ) )  =  1 >  0E  6 
AO  CONTINUE 

BEGINING  OF  THE  SEARCH  FOR  BEST  SUBSET 

DO  99  M=MS  t ML 
NR=N -M 

DO  70  1=1 > NR 
IND( I )=1 
INDB< I ) =1 
70  CONTINUE 

DO  80  J=1 f M 
IND< J+NR) =0 
INDEK  JFNR>=0 
80  CONTINUE 
IE=~1 
MF’=NR+1 
BR=AE  <  MP  ) 

KP=1 

KPT=KF*TH 

MI*M+1 

ETAN-l .OE  6 

DO  90  I T  =  1  f  M I 

IF<IT*EG»»1)60  TO  82 

KP=(KP*(NRfIT  -2) )/< IT-1  ) 

KPT  -  KPTTKP 

82  DO  86  K=1 » KP 

IF  < IT » EQ » 1 )G0  TO  83 
CALL  COMBO < IND f N f NR f NF ) 

83  CALL  INTG<AfCOLfYf INDfIPfNf NSfLFI f IE) 
CALL  AHERI V  < A f  COL  f  SOL  fYfETfNSf KN ) 
ETA=CDABS<ET) 

IF(ETA,GT.ETAN)G0  TO  86 

ETAN=LTA 
DO  84  1=1 fN 
INDB<I)=IND<I) 

84  CONTINUE 

DO  85  J= 1 f  M 
SOLI < J)=SOL( J) 

85  CONTINUE 

86  CONTINUE 

87  IF  <  ETAN .  I  E . BR )G0  TO  92 
80  MP=MP  f 1 

BR=AE ( MP ) 

90  CONTINUE 
92  J=0 

DO  93  L=1 »N 


***arf**— •• —  . .  •. ^ 


IF<INDB<L>  .EQ.DGO  TO  93 
J=Jfl 

JP< J)=IP(L) 

93  CONTINUE 

DO  95  LP=lfJ 
L JP  <  LP ) = JP ( LP ) 

JJP  <  LP ) =LP 
DO  94  KF’-l » J 

IF ( L JP < LP ) . LE . JP ( KP ) ) GO  TO  94 
L JP  <  LP ) = JP  <  KP ) 

JJP  <  LP ) =KP 

94  CONTINUE 

JP< JJP<LP> >=1E  6 

95  CONTINUE 

APPLYING  WINDOWS  AROUND  SELECTED 
BASIS  FUNCTIONS 

NO=NO+l 

ni=n.t/2 

D2=ni 

I F  ( NO  » GE . 2 ) GO  TO  505 
D3=D  <  L  JP  <  1 ) )  -D 1 
D4=D<LJP<2> >—Dl 
GO  TO  506 
505  MS--2 
ML=2 

IF ( DABS < D ( L JP < 1 ) ) -D < L JP < 2 ) ) ) .LE. .002D0)G0  TO  550 
D3--D<LJP<1 )  )-2*Dl 
D4=D(LJP<2) >-2*Dl 
50A  IF ( D3 . GE . O . DO ) GO  TO  520 
D3=0 . DO 

IF (D4~D3>525>525»520 
525  D4-D4f2*IU 

520  IF <D4~D3~4*ni >501 .502,503 

501  D4=D4i3#Di 
GO  TO  503 

502  D4«D4+D1 

503  IF (N0.GE.7JG0  TO  93 

IF  <  ETAN . GT . 1 . 0E--10  )  GO  TO  3 

98  CONTINUE 

99  CONTINUE 

550  EDI < I J>=D(LJP( 1 ) ) 

ED2< I J)=D<LJP(2) ) 

WRITE (6, 300) FBI < I J) ,FD2( I J) 

500  CONTINUE 


EVALUATION  OF  MEAN  AND  VARIANCE  OF 
THE  ESTIMATES 


AD-A09Q  126  RHODE  ISLAND  UN IV  KINGSTON  DEPT  OF  ELECTRICAL  ENGIN— ETC  F/G  9/2 

APPLICATIONS  OF  efficient  subset  selection  TO  digital  filtering — ETC<U) 
AUG  80  J  KHORAHMI *  D  TUFTS  AFOSR-77-3174 

UNCLASSIFIED  AFOSR-TR-80-095A  NL 
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CALL  MARIAN ( EDI »ED2 » DMEAN1 > DMEAN2 > MARI » MAR2  > NRUN > 
WRITE(A>290) 

DO  510  111=1 >NRUN 

Enil<IIl)=DABS<EBl<IIi)-Bl) 

En22<IIl)=nABS(En2<IIl)-B2) 

WRITE  <  A>  300)ED1  <II1)>ED11(II1)  >ED2(  1 1 1 )  > EB22(  III) 

510  CONTINUE 

BMER1=BABS(B1-DMEAN1 ) 

DMER2=DABS  <  B2-BME AN2 ) 

STDl=nSQRT (MARI ) 

STD2=DS0RT<UAR2) 

WRITE(A>310)DMEAN1>DMEAN2 
WRITE <  A  > 330 ) BMER 1 > DMER2 
WR I TE  <  A  > 320 ) MAR 1 > V AR2 
WRITE <A> 340 )STD1>STD2 
100  FORMAT <11 >  3I2>  2F10 . 0 ) 

105  FORMAT <213*212) 

110  FORMAT < I3>4F10.0) 

120  FORMAT <5X> 'D('>I2>')='>E14.A> 

130  FORMAT </> 5X > 'C ('» 12 >' )  =  <>E12.4>'  >  '  * 

*E12>4>'  )') 

140  FORMAT <//> 5X> 'RSS  FOR  FULL  BASIS  ='>E12.4) 

150  FORMAT < // > 5X > ' N= ' > 12  » '  M='>12) 

1 AO  FORMAT <  /  > 5X » 'E< ' >I2> ' )=' >E12.4) 

170  FORMAT <  /  > 5X > '  E  < ' >  12  > ' )  = ' >  E12 . 4  > '  NO.  '>12 
*>'  SMALLEST  ERROR') 

180  FORMAT  < / >  5X  >  'C<  '  >  12  >  '  )  =  (  '>E12.4>  '  >  '>E12.4>'  )') 

190  FORMAT  < // >  5X  > ' ERROR= '  >  E 1 2 . 4  ) 

200  FORMAT <//>5X> 'NO.  OF  SUBSETS  EXAMINED  =  '>14) 

210  FORMAT  < // > 5X  > ' N= ' > 12 ) 

220  FORMAT <2I2>2F10»0) 

230  FORMAT ( / >  5X > ' KN= ' > 1 3  >  3X  > ' NRUN= '>13) 

280  FORMAT </>5X> 'SIR  =  ' >E10 ♦ 4  > 3X  > ' SNR  •=' >E10.4>3X 
*  > ' NM AR= ' , E 1 0 . 4 ) 

290  FORMAT ( // >  9X  > ' FREQ1 ' >7X> ' FREQ  - ERR0R1 ' >7X> ' FREQ2 ' 

*>7X> ' FREQ-ERR0R2 ' ) 

300  FORMAT  < / >2X  >4<3X>E12.A) ) 

310  FORMAT < // > 5X > 'MEAN1  =  ' >E12. A>7X> 'MEAN2  =  '>E12.A) 

320  FORMAT ( / > 5X > ' MAR I ANCE 1='>E12.A>5X> ' MARIANCE2= ' >E12 . A > 
330  FORMAT </>5X> 'MERR0R1= ' > El 2 . A > 7X > ' MERR0R2= ' > E 1 2  *  A ) 

340  FORMAT (/> 5X >' STDIM1  =  ' > E12 . A > 5X > ' STDIM2  =  '>E12.A> 
350  FORMAT </>5X> 'Al  =  ' >E10.4>3X> ' A2= ' > El 0 . 4 ) 

3A0  FORMAT ( / > 5X > 'Fl= ' > E10 . 4 > 3X > ' F2= ' >E10 . 4 > 

STOP 

END 


SUBROUTINE  INTG  IS  FOR  EMALUATION  OF  ELEMENTS  OF 
A  AND  B  MATRICES  FOR  ACM  SUBSET  TO  BE  EXAMINEDf 
EQUATIONS  (4.2.14)  THRU  <4.2. 1A) 
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SUBROUTINE  INTG  <  An  COL » Y » INDf  IP»N»NS f  LFI » IE ) 
COMMON  PI,A1»A2»B1 f B2»D1 fD2»D3»D4»N1 »N2»KN 
C0MPLEX*16  A<  40f40) fC0L<40>  f Y< 100) 

REAL*8  D(40) f AD<  40) fX  f  Y1 fY2fY3fA1 f A2f  B1 f  B2  fPI 

REALMS  D1fD2fD3fD4 

INTEGER  INK (40) f IP<40) » JP(40) 

IF  <LFI ♦ EQ » 1 )G0  TO  2 
D<1)=D1  ^ 

BO  1  L=2,N  \  ’ 

1  D(L)~B(L~1 )+Dl 

GO  TO  5  7"  " 

2  D  <  1 )  --D3 
D ( 6 ) =D4 

DO  3  L=2 f 5 
D<L>=D<L-1>+D1 

3  D<L+5)=D<L+4 )+Dl 

5  IF  < IE . EQ  »  0 )G0  TO  30 
IF ( IE . EQ .  -1 )G0  TO  25 
NS=N  -1 
DO  20  L=1 f N 
IF<L .NE ♦ IE )G0  TO  20 
DO  15  M-LfNS 
15  D<M)~D<M+1) 

20  CONTINUE 
GO  TO  30 

25  K”0 

Kl  =  l 

DO  26  K2=l t N 
IF<IND(K2)  .EQ.DGO  TO  26 
AD  <  K 1 ) =D  < IP  <  K2 ) ) 

K=K+1 

K1=K1+1 

26  CONTINUE 

DO  20  K3~l fK 
D(K3)=AD(K3) 

JP<K3)=K3 
DO  27  K4-1 »K 

IF ( D(K3 ) *LE . AD< N4 ) ) GO  TO  27 
D(K3)=AD<K4) 

JP<K3)=K4 

27  CONTINUE 

AD< JP<K3) >=1 .OE  6 

28  CONTINUE 
NS=K 

30  DO  60  I  =•  1  f  NS 
DO  60  J=I fNS 
XSS<D<J)~D(I>  )#2.0D0*PI 
IF  <X)40f50f40 
40  CONTINUE 
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A  < I >  J  >  =  < 1 . OnO-CDEXP <  DCHPLX ( 0 . DO » X*KN ) > ) 

A  <  I » J  )=A  (Ifj)/<1»  DO~CEiEXP<  DCMPLX  <  0 » DO » X )  ) ) 

A  <  J  »  I  >  =DCON JG <  A  <  I  »  J  > > 

GO  TO  60 

50  A ( I » J  >  =  1 • ODO*KN+  » 00001 
60  CONTINUE 

no  80  K=1,NS 
coL<K)=(o.no»o.no> 

Y3=-2.no*Pi*n<K> 
no  80  J=1»KN 

CQL<K)=COL<K)+Y< J>*CDEXP<BCMPLX<0.n0,Y3*< J~1 >  >  ) 
80  CONTINUE 
RETURN 
ENO 


SUBROUTINE  AHERIV  IS  FOR  SOLOING  THE  LINEAR 
EQUATIONS  AC=B  AND  CORRESPONDING  ERROR  FUNCTION 

SUBROUTINE  AHERIV < A » COL , SOL » Y » ET  *K » KN  > 

COMMON  PI  * A1 » A2»  B1 » B2 1 M  » D2 » D3»  D4 » N1 » N2 
REAL*8  A1 *  A2 » B1 » B2 , PI » C » ON ( 50 ) > Cl t C2 
COMPLEX* 1 A  A ( 40 1>40 ) » COL  <40 ) »RL  < 40r  40 ) t SOL ( 40 ) 
COMPLEX* 1 6  XXI  *  Y  < 1 00  >  t E 1  *  E2  * E3  > ET t XX 
DO  5  1=1, K 
DO  5  J=.t  ,  K 

5  RL  <  I » J )  =DCMPLX  <  0  .  ODO » 0 .  OHO ) 

no  60  1=1 fK 

IP1=I+1 
I M 1  =  I  - 1 

1F<I.NE.1)G0  TO  30 
DO  10  J-- 1  »  K 
10  RL  <  1 ,  ,J ) =A  < 1 » J  > 

DO  20  J=2»K 

20  RL< Jfl >=DC0NJG(RL<1 , J> )/RL(l , 1 > 

GO  TO  60 
30  BO  AO  J=I*K 

RL  <  I  ,  J  >  =  A  ( I  ,  J  > 
no  AO  M=1 ,  IM1 

AO  RL(  I ,  J>=Rl.  <I,J>~RL<I,M)*RL(M,J> 

IF ( I ♦ EQ  »K ) GO  TO  60 
DO  45  J=IP1,K 

45  RL  <  J  *  I ) =DCON JG <  RL  < I ,  J  >  ) /RL ( I » I ) 

60  CONTINUE 
80  DO  100  1=1 ,K 
IM1=I- 1 
SOL  < I )=COL  < I ) 

IF  < I . EQ » 1 ) GO  TO  100 
no  ?0  L=1 » IM1 

90  SOL  ( I )  =SOL  ( I  >  -SOL  ( L  )  #RL  <  I  r  L  > 
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100  CONTINUE 

DO  140  1=1 » K 
J-K-I+l 

IF < J ♦ EQ »K )G0  TO  130 
JP1=J+1 

no  120  L= JP1 » K 

120  SOL  <  J )  =S0L  <  J  >  -SOL  <  L  )  *RL  ( J » L  > 

130  SOL  <  J ) =S0L ( J ) /RL  <  J ,  J ) 

140  CONTINUE 

EVALUATION  OF  ERROR,  EQUATION  (4.2.19) 


E2=0.0ti0 
Ei=o.ono 
no  15  11=1, K 

E2=E2+S0L  (ID  *DCON  JG  <  COL  <I1>) 

15  CONTINUE 

HO  16  13=1, KN 
El=El+CriABS(Y<I3)  >**2 

16  CONTINUE 
ET-E1-E2 
RETURN 
END 


SUBROUTINE  COMBO  IS  FOR  GENERATING  STRING 
OF  ONES  ANH  ZEROS  FOR  PROPER  BELETION  OF 
BASIS  FUNCTIONS 

SUBROUTINE  COMBO ( JN , N , M , NF ) 

DIMENSION  JN ( 40 ) 

IONES=M 

NF=0 

1  NP=1 

2  CONTINUE 

IF< JN(NP) ,EQ.O)GO  TO  10 
IF  <  NP . EQ . N  >G0  TO  101 
JN  ( NF' )  =0 
NP=NP*1 
I0NES=I0NES~1 
GO  TO  2 
10  JN<NP  >  =  1 

I0NES=I0NES+1 
IF < I ONES . EQ »M)GO  TO  100 
GO  TO  1 
101  NF=1 
100  RETURN 
END 


•  0 
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SUBROUTINE  GAUSS  AND  RANDU  ARE  FOR 
GENERATING  GAUSSIAN  NOISE 

SUBROUTINE  GAUSS< IX»S,AM,V> 

A=0.0 

DO  50  1=1 » 12 
CALL  RANDU < IX»IY»Y) 

IX=I  Y 
50  A=A+Y 

V=<  A~6  »  0  )#S+AN 

RETURN 

END 

SUBROUTINE  RANDU < IX , IY » YFL > 

IY=IX#65539 

IF(IY)5»6*6 

5  IY=IY+2147-4836-47+l 

6  YFL=I Y 

YFL=YFL* .  -465661 3E-9 

RETURN 

END 


SUBROUTINE  VAR I AN  IS  FOR  EVALUATING  MEAN 
AND  VARIANCE  OF  THE  FREQUENCY  ESTIMATES 

SUBROUTINE  VAR I AN <  EDI » ED2 , DMEAN1 , DMEAN2 , VAR1 r VAR2 > NRUN ) 

REAL*8  EDI < 100) » ED2 ( 100  >  r  DMEAN1 »  DMEAN2  r VAR1 

REAL*8  VAR2 > XI » X2 

VAR 1=0. DO 

VAR2=0 . DO 

DMEAN1=0 . DO 

DMEAN2=0 . DO 

DO  10  1  =  1  » NRUN 

DMEAN1=DMEAN1 +ED1 < I ) 

DMEAN2  =DMEAN2+ED2 < I ) 

10  CONTINUE 

DMEAN1 =DME AN1 /NRUN 
DMEAN2=DMEAN2/NRUN 
DO  20  J=1 » NRUN 

VAR1=VAR1+<ED1 < J)-DMEAN1 >**2 
VAR2=VAR2+ ( ED2 ( J ) -DMEAN2 ) **2 
20  CONTINUE 

VAR 1=VAR1 /NRUN 

VAR2=VAR2/NRUN 

RETURN 

END 


SUBROUTINE  BELT  IS  FOR  DELETING  PROPER 
ROW  AND  COLUMN  OF  A  AND  B  MATRICES  FOR 


EVALUATION  OF  El  ERROR  FUNCTIONS 


SUBROUTINE  BELT  < A » COL * S *  DCOL *  N  *  NS » IE  > 
C0MPLEX*1&  A<40»40) » COL (40) *S<40*40) »DCOL<40> 
IEM=IE-1 

IF ( IE  »EQ . 1 )G0  TO  30 
DO  10  1*1* IEM 
DCOL  < I >  =COL  < I ) 

DO  10  J=I  * IEH 
S<I*J)=A<I* J> 

S< J*I>=A< J*I> 

10  CONTINUE 

DO  20  1=1* IEM 
DO  20  J=IE »NS 
S< I  * J)=A< I » J+l > 

S  <  J » I ) =A<  J+l *  I > 

20  CONTINUE 

IF<  IE . EQ.N)GO  TO  50 
30  DO  40  I  =  IE *NS 

DCOL ( I )=COL ( I +1 ) 

DO  40  J=I »NS 
S<I*J)=A<I+1»  J+l) 

S( J.I)=A<J+1*I+1) 

40  CONTINUE 
50  RETURN 
END 
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